Direct path ranging experiment

Our inability to resolve sub-horizontal tectonic displacement on the seafloor is addressed by a path-ranging method based on acoustic travel-time measurements16,20,21,22,25. Distance changes between seafloor instruments are estimated from two-way travel times and sound speed of water. Sound speed along the entire ray-path is approximated by the geometric mean of the sound velocities measured at both endpoints. In October 2014, an offshore geodetic network of intercommunicating transponders was installed in 800 m water depth where the NAF trace is identified in high-resolution multibeam seafloor bathymetry maps26,27 acquired by autonomous underwater vehicles (Fig. 2) and in 3.5 kHz seismic profiles28. The deployment site was selected based on the existence of a linear scarp of the NAF clearly visible in the bathymetry as a proxy for maximised strain release along the NAF29,30. Ten acoustic transponders, four from the Ocean Geosciences Laboratory in Brest, France (station names starting with F) and 6 from GEOMAR’s GeoSEA array in Kiel, Germany (station names starting with G), were installed25 in October 2014 and remained fully operational until May 5, 2017. Each set of stations, manufactured by Sonardyne Ltd (UK), used a different center-frequency (F: 22.5 kHz, G: 17 kHz), so the F- and G- stations communicated only with F- and G-stations, respectively. Pairs of stations, <100 m apart, shared common baselines and all stations monitored the temperature, pressure, and tilt (Fig. 3a, b, Supplementary Figs. 1–3) at the transponder site along with the two-way travel time between them (Fig. 3c).

Fig. 2 Map of the geodetic network on the seafloor. Baseline shortening or lengthening rates, based on the total deformation measured at the end of the deployment, are colour coded. A constant strain of 4.5 × 10−6 was subtracted to each baseline in order to correct for the inferred constant salinity decrease of 0.002 PSU yr−1. Transponder locations are shown as black circles. F- and G- stations only communicate within their respective network. High-resolution shaded bathymetry map26,27 is shown with a light source from the South Full size image

Fig. 3 Measured parameters and estimated baseline for station pair G2-G5. Black lines indicate monthly medians. a Temperature. b Pressure. c Acoustic one-way travel times between the transponders for the forward- and backward measurement. d Estimated sound speed from temperature and pressure (using a constant salinity). e Resulting baselines for both measurement directions. Orange lines indicate the monthly mean baseline lengths with the inferred linear salinity decrease rate of 0.002 PSU yr−1 Full size image

Baseline estimation

We use 650.000 two-way travel time measurements based on two soundings, respectively. Each station either interrogates the other stations of its kind or acts as a replying instrument, thus the experiment forms an autonomous intercommunicating network on the seafloor, rather than observations at individual positions. The acoustic distance is then calculated from the average sound-speed in water multiplied by the one-way travel time.

We estimate the sound-speed in water from temperature and pressure measurements32, assuming a constant salinity, similar to other path-ranging experiments20,33,34. Figure 3 shows the time series of the measured parameters together with the calculated sound speed and baselines for the southwest-northeast fault crossing baseline G2-G5.

For the final strain estimates (Fig. 4), we subtracted a mean strain of 4.5 × 10−6 estimated from all baselines not crossing the fault (corresponding to a baseline decrease of 4.5 mm for a 1000 m long baseline during 2.5 years). Subtracting a constant strain value, instead of a constant increase or decrease in the baseline lengths, accounts for the fact that linear trends in sensor drift or water parameters are linearly related to strain. In addition, considering strain makes the observation independent from the baseline lengths. The baselines show a consistent behaviour (Fig. 2, Supplementary Figs. 4, 5 and Supplementary Table 1), with a long-term resolution better than 8 mm for all baselines (Fig. 3), and strain is smaller than 1 × 10−5 for all baselines (e.g. less than 1 cm deformation on a 1 km-long measurement distance). For instantaneous baseline changes, the resolution capability is ~2 mm, in particular, if an event would influence different baselines. We interpret the nominal lengthening of all baselines as a result of subtle changes in the physical properties of water, such as a salinity decrease of 0.005 practical salinity units (PSU) during the 2.5 years of deployment (0.002 PSU a−1). Alternatively, one could consider a total drift in all pressure sensors of −0.242 kPa or a total drift in all temperature sensors of −0.00126 °C. In particular, a pressure drift of −0.242 kPa during the deployment cannot be distinguished from a local uplift of 2.4 cm. Because the Sea of Marmara is an extensional step-over or pull-apart structure system including substantial subsidence35 with the basement imaged at a depth of ~4.5 km below the Central High and therefore close to the geodetic network, uplift is unlikely36. The data suggest an absence of vertical movement and the remaining small baseline changes originate from measurement uncertainties of the pressure and temperature sensors. The subtle changes of water parameters such as the temperature increase of 0.002 °C on March 2016 (Fig. 3a) might be a sensor artefact since they are not compensated by pressure or travel time resulting in an apparent ~3 mm baseline length change (Fig. 3e). However, since baseline estimates are based on the equation distance = time × velocity, travel times and water sound-speeds must be jointly accurately known. This problem is similar to the hypocenter-depth-velocity dependency in earthquake location techniques. This is the reason why the network was designed to measure a high number of baselines across the fault to allow isolating effects of sensor drift from baseline changes.

Fig. 4 Strain of all baselines. Strain (monthly medians) estimated from the measured travel times and calculated sound speeds based on pressure, temperature and a linear salinity decrease of 0.002 PSU yr−1. Baselines not crossing the fault (e.g. subparallel to the NAF) are shown in black; for a simple, east-west oriented, right-lateral strike-slip fault, baselines oriented southwest-northeast (red lines) should lengthen and those oriented northwest-southeast (blue lines) should shorten. Predicted strain rates correspond to a strike-slip movement reaching the seafloor Full size image

Estimation of slip on the fault from baseline data

We compared the observed baseline changes with a vertical west-east trending strike-slip model crossing the network. We used a least square inversion to determine the slip rate of the fault which minimizes the differences between the observations and the strike-slip fault model25. This approach implicitly includes the assumption that baselines located on one side of the fault (i.e. not crossing the fault) are not changing. From baselines crossing the fault we found an optimal rate for strike-slip movement of 0.80 ± 1.25 mm yr−1. This suggests that the surface fault slip rate across the network is close to zero and consistent with the results from the first six months of deployment25. Analysis of the geodetic data during the first six months of the deployment resulted in an upper bound on the slip rate of 6 mm yr−1 only25. The seven-fold increase of fault slip resolution clearly demonstrates the need for long-term deployments in order to resolve tectonic processes.

Modelling of strain and locking depth

Next, we model strain for vertical strike-slip faulting. We use an analytical half-space solution, based on the elastic dislocation theory for an infinitely long vertical strike-slip fault37. Since the geodetic network is located on low rigidity sediments, we model the dependency of strain in the presence of a low rigidity layer (=weak layer). We use a simple model consisting of two horizontal layers with a rigidity contrast37 to investigate the dependency of strain on creep below given depths, rigidity and varying thickness of the overlying sedimentary layer. The slip rates estimated for the NAF from onshore observations range between 15 and 27 mm yr−1 10,38,39 from GPS observations and between 15 and 19.7 mm yr−1 from mass deposit considerations40,41. We model fault creep with a rate of 20 mm yr−1 (in-between the GPS and geological estimates) for different faulting depths, above which the fault is locked and below which it creeps. The modelled strain is linearly dependent on the inferred slip rate of 20 mm yr−1 37. We model scenarios for creep below 3 km (corresponding to the thickness of the basin sediments36) and hence strain in the pre-kinematic basement rocks (Fig. 5a) and for creep below 4.5 km depth (corresponding to the depth of the crystalline basement) (Fig. 5b). The depths of the geological units are known from a seismic profile passing in ~5 km distance south-east of the geodetic network36. The results show that the overlaying low rigidity sedimentary layers focus the strain close to the fault (Fig. 5a, b). Fault creep below 3 km (Fig. 5a), corresponding to slip in the pre-kinematic basement rocks36, is clearly above the strain rate sensitivity of the geodetic network. The strain rate sensitivity of the geodetic network is 1.6 × 10−6 yr−1 corresponding to the inverted strike-slip movement of 0.8 mm yr−1 considering the 500 m coverage of the geodetic network perpendicular to the fault (Fig. 5). Modeling the minimal fault slip inferred from onshore geodetic observations (16 mm yr−1) results in 20% less strain and would still have been above to the sensitivity of the offshore geodetic network, in particular, due to the existence of weak shallow layers (Fig. 5a). From the modelling we find the locking depth of 3 km as the most conservative estimate. The strain rate induced from 20 mm yr−1 creep in the crystalline basement (Fig. 5b), located below 4.5 km36 depth, results in a signal exceeding the strain rate sensitivity for sedimentary layers thicker than 1 km.

Fig. 5 Modelled strain rates for a vertical strike-slip fault. The strain for homogenous and layered half space, based on elastic dislocation theory, was estimated using the analytical solution for a horizontally layered half-space and a vertical strike-slip fault37. a Model for creeping below 3 km depth corresponding to slip in the pre-kinematic basement rocks36 and below. Above the locking depth, the fault is fully locked. Below the locking depth, the fault creeps at 20 mm yr−1. The ratio of rigidity between the lower and upper layer is 4. The strain rate sensitivity of the geodetic network is 1.6*10−6 yr−1, corresponding to strike-slip movement of 0.8 mm yr−1 considering a 500 m extension of the geodetic network perpendicular to the fault. b Model for slip below 4.5 km in the basement rocks36. c Creep below 16 km depth corresponding to interseismic deformation. d Dependency of rigidity contrast on the strain rate for a 4.5 km weak layer and creep below 5.5 km. Rigidity contrast from empirical relations is in the range between 8 and 16 (Supplementary Table 2). e Sketch showing the modelling setup and parameters. μ 1 and μ 2 are the rigidities of the shallow layer and below Full size image

Slip at depths below 16 km (Fig. 5c), corresponding to the inter-seismic deformation, results in a small strain signal and is still in line with the measurement uncertainty of strike-slip faulting of 0.80 ± 1.25 mm yr−1. In the last step, we modelled the dependency of strain on rigidity contrasts between the upper and lower layer (Fig. 5d). We estimated the rigidity from empirical relations from a seismic P-velocity profile36 (Supplementary Table 2). Despite the uncertainties of the empirical relations we find rigidity ratios clearly exceeding 10 for the material above and below the basement, so the calculated strain rate of strike-slip faulting on the seafloor is little dependent on these large values (Fig. 5d). From the strain rate estimates, we conclude that the strain rate sensitivity of the geodetic network is likely sufficient to resolve fault movement above the basement (Fig. 5a, b) and the model shows that the maximal resolution is reached for slip occurring at depths below 5.5 km (Fig. 5d).

With the possibility of distributed strike-slip across a few kilometre-wide zone of faults at the seafloor26,30,42,43, we might not have captured the complete possible slip. However, deformation is clearly focussed beneath the geodetic network since the fault trace can be unambiguously identified in the bathymetry26,30. In particular, there is a clear 3.5 km right-lateral offset of a ridge between the baselines, indicating the location of the geodetic network above the main zone of surface deformation of the NAF29.

Local seismicity

To better detect small-magnitude events indicative of a creeping behaviour, two small aperture OBS arrays were deployed in the vicinity of the geodetic stations and close to the NAF: a 5 km wide array during five months and a 12-km wide array for the next 12 months (Figs. 1 and 6). Such small aperture OBS arrays are significantly more sensitive to low-magnitude (from 0 and up) and shallow seismicity than larger aperture OBS arrays (e.g., 10 km station spacing) which have typically a magnitude of completeness of 113. The OBS detections were complemented with phase picks from the land stations of the Kandilli Observatory and Earthquake Research Institute (KOERI). In 17 months, only 45 events with local magnitudes between 1.4 and 4.2 were detected and none are located near the geodetic experiment (Fig. 6). These events are 9 to 25 km deep and slightly offset (~2 km) from the NAF surface trace. The offset can be explained with the increased strain and hence stress accumulation in the vicinity of the locked fault37 (Fig. 5c) and gas migrations at depth42,44. The local observation is in-line with the sparse seismicity recorded by onshore stations during the last decade along this segment of the NAF (Fig. 1) and, as the geodetic observation, also concurs with its quasi-locked status.