# Creep-dilatancy development at a transform plate boundary

Nabil Sultan, Shane Murphy, Vincent Riboulot, Louis Géli

► To cite this version:

Nabil Sultan, Shane Murphy, Vincent Riboulot, Louis Géli. Creep-dilatancy development at a transform plate boundary. Nature Communications, 2022, 13, 10.1038/s41467-022-29558-0. insu-03685390

HAL Id: insu-03685390

https://insu.hal.science/insu-03685390

Submitted on 2 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Creep-dilatancy development at a transform plate boundary

Nabil Sultan<fnref n="1" />, Shane Murphy<fnref n="1" />, Vincent Riboulot<fnref n="1" /> & Louis Géli<fnref n="1" />

How tectonic plates slip slowly and episodically along their boundaries, is a major, open question in earthquake science. Here, we use offshore in-situ sediment pore-pressure acquired in the proximity of the active offshore Main Marmara Fault and onshore geodetic time-series data set from a single GPS station to demonstrate the pore-pressure/deformation coupling during a 10-month slow-slip event. We show that pore pressure fluctuations are the expression of hydro-mechanical process affecting the deep seismogenic zone and indicate that small disturbances in geodetic data may have important meaning in terms of transient deformations. These results have major implications in understanding the spatial impact of slow-slip processes and their role in earthquake cycles. We demonstrate that piezometers measuring along a transform fault can help define the time scale regulating the coupling between slow-slip events and earthquake nucleation process.

<footnote n="1">¹Geo-Ocean UMR6538, Ifremer, CNRS, UBO, UBS, 29280 Plouzané, France. ☎email: nabil.sultan@ifremer.fr</footnote>
That tectonic plates may slip slowly and episodically along segments of their boundaries is one of the most intriguing, unexplained observations in solid earth geophysics, as it reveals a continuum of transient deformation along active faults, ranging from seismic rupture to aseismic events¹,². While the major role and influence of pore pressure has long been proposed³⁻⁵, it has never been observed because of the difficulty in linking deep processes at seismogenic depths with near-surface pore-pressure measurements. Slow slip events (SSEs) have mostly been studied using on-shore Global Positioning System (GPS) data⁶, where the discovery of fluid flow based on direct, in situ observations from the seafloor along the subducting plate interfaces have led to the hypothesis of a causal relationship between SSEs and changes in fluid activities at the fault zone⁷⁻⁹.

The Sea of Marmara—south of Istanbul—has become, after 20 years of investigations, one of the most well-known submarine domains on earth. Studies have now largely converged to a scheme with a master fault—the Main Marmara Fault¹⁰ (MMF; red line in Fig. 1A)—taking 75% of the right-lateral geodetic slip rate (~20 mm/year) between the Eurasian Plate and the Anatolian block¹⁰⁻¹³. The western segments of the MMF are characterized by a higher level of background seismicity compared to the eastern segments. The common occurrence of small repeating

earthquakes¹⁴ suggests that the MMF has a deep creep component in the western Sea of Marmara, but not in the central nor in the eastern segments¹⁵. Based on submarine, acoustic ranging data, aseismic creep rates of 9–16 mm/year¹⁶ were measured in the Western High area (WH in Fig. 1A), where active seeps¹⁷⁻¹⁹ were found in association with an active mud volcano (MV) located <1 km away from the MMF trace (Fig. 1B, C). The largest gas emissions along the MMF were found to occur within a 900 m to 1000 m zone on either side of the MMF²⁰, which encompassed the MV and have been interpreted as being the trace of a high-permeability damage zone²¹.

Here we analyze offshore in situ sediment pore-pressure acquired from the active MMF and onshore GPS time-series data. Our results indicate that a pore-pressure/deformation coupling exists during a 10-month slow-slip event. These observations demonstrate the role of dilatancy in regulating tectonic slip, reactivating an adjacent MV and far-field tectonic deformations.

## Results

In this context, two differential piezometers²² measuring pore pressures in excess of hydrostatic pressures (“Methods”) were deployed from Oct 2013 to Nov 2014²³ within the MMF valley

<figure><img src="image_3.png" /><figcaption>**Fig. 1 Study area.** A General map showing the study area, the onshore GPS stations (dots indicating from west to east: CANA, YENC, TEKR, BALK, and ISTN stations), and the 2011EQ, 2013EQ, and 2014EQ earthquake epicenters (black stars). The thick black line¹¹ indicates the trace of the northern strand of the North-Anatolian Fault. Background map boundaries are provided by ArcGIS Data and Maps. Redistribution rights are granted. https://www.esri.com/content/dam/esrisites/en-us/media/legal/redistribution-rights/redist-rights-2021.pdf. B Zoom of bathymetry showing piezometers PZN and PZS and their positions with respect to the Main Marmara Fault (MMF) and its damaged zone (Supplementary Fig. 1), the mud volcanoes (MVs), and gas seep areas¹⁷. C 3D diagram including interpreted seismic lines extracted from 3D HR seismic data and indicating the position of PZN and PZS.</figcaption></figure>
and close to the MV, respectively (PZS and PZN in Fig. 1), to
investigate the coupling between the hydrogeological and tectonic
systems along the WH segment of the MMF. The PZS and PZN
water depths are 667 and 672 m, respectively.

During the monitoring period, two notable earthquakes
occurred at epicentral distances from PZN of 12.5²⁴ and
209 km²⁵, respectively: the 11 km deep, Mw 4.8 27/11/2013
(2013EQ) nucleated along the MMF in the area where repeating
earthquakes have been documented, while the 9 km deep, Mw 6.9
24/5/2014 earthquake (2014EQ) struck the Ganos segment of the
NAF in the Thracian Sea (epicenters in Fig. 1A).

In parallel, we analyze 7-year period GPS data²⁶ using publicly
available data from Nevada Geodetic Laboratory-NGL (http://
geodesy.unr.edu/index.php). GPS time series data are provided in
the IGS14²⁷ reference frame and in the plate fixed reference frame
by considering the plate on which the station lies. For this last
case, the plate motion trends are removed for the horizontal
component time series. The accuracy of an observed position in
the north-south, east-west, and in up-down directions depends
significantly on the duration of the observing period and on the
accuracy of the GPS satellite orbits used during processing²⁸. In a
previous study, an accuracy analysis of relative positions of per-
manent GPS stations in the Marmara Region has shown that the
root mean square error is within 1 mm for the north-south,
east-west components while it is between 2 and 3 mm for the
up-down direction¹⁶. For the present study, we used the STL
(seasonal-trend decomposition based on LOESS)²⁹ procedure in
R³⁰, which is a filtering procedure that decomposes a time series
into three components: trend, seasonal, and remainder parts.
Figure 2 shows a decomposition plot of TEKR geodetic data with
a trend of the northern displacement (with respect to fixed EU
plate) strongly affected by the 2013EQ. The vertical displacement
at TEKR also indicates a perturbation of the trend between
2013EQ and 2014EQ.

The horizontal displacement (north with respect to fixed EU
plate)⁶, recorded by the four westernmost stations, indicates a

general evolution to the south in agreement with the geodetic strain
rates within the region³¹ (Fig. 3A). The non-corrected vertical (up)
GPS displacement indicates seasonal variations (Fig. 3A). Super-
imposed on this general trend, displacements (north and up) at
station TEKR seem disturbed (arrows in Figs. 2 and 3) during a
period following the 2013EQ. The eastern horizontal GPS dis-
placement from the Anatolian plate evolves smoothly to the west
independently of 2013EQ (Supplementary Fig. 2).

Excess pore-pressure data (i.e., above hydrostatic) measured at
7.08 mbsf (PZN-P6) and 6.28 mbsf (PZS-P5) show two opposite
trends (Fig. 3B and Supplementary Figs. 3–6). The PZN-P6
pressure increases over a period greater than 3.5 months starting
around 1 week after the 2013EQ. Pore-pressure and temperature
data from PZN-P6 and PZN-P5 suggest a diffusion–advection
process controlling the thermal and pressure gradients (Supple-
mentary Figs. 7, 8, and 10). During this period, the pore pressure
repeatedly reaches hydro-fracture conditions (Supplementary
Fig. 11). After this, a smooth decay of the pressure is followed by
a 3-month period of negative pore pressure (gray dashed area in
Fig. 3B). The sensor positioned 0.8 m above P6 was the only one
to record significant pressure changes during the monitoring
period (Supplementary Fig. 3). At PZS-P5 station, negative pore
pressure develops 1 week after the 2013EQ and lasts for
2.5 months (Fig. 3B). The P5 sensor was the only one to record
significant pressure changes during deployment (Supplementary
Fig. 5). The zoom diagrams in Fig. 3C show in detail the pore
pressures generated by 2013EQ and 2014EQ during the mon-
itoring period. The five other sensors of PZN and PZS were also
disturbed by those two events (Supplementary Figs. 4–6).

**Near-field pore-pressure data versus far-field geodetic dis-**
**placements.** Comparison in Fig. 4A between the Tekirdag daily**
**mean temperatures (from climexp.knmi.nl) and the vertical dis-**
**placement from TEKR station confirms the seasonal origin of the**
**geodetic vertical component, where the redistributions of con-**
**tinental water mass deform elastically the Earth crust similar to**

<figure><img src="image_4.png" /><figcaption>**Fig. 2** Decomposition plot of TEKR GPS geodetic data. **A** Northern displacement and **B** vertical displacement. The unit on the vertical scale is meters. The horizontal scale uses the convention of day, month, and year. The timing of the three earthquakes (2011EQ, 2013EQ, and 2014EQ) are displayed by dashed vertical lines.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Time series data. **A** Displacement time series (from left to right: north and up components corrected for seasonal fluctuations and non-corrected up component) for five sites surrounding the Sea of Marmara. From west to east: CANA, YENC, TEKR, BALK, and ISTN stations (locations in Fig. 1A). The three vertical dashed lines indicate the timing of three neighboring earthquakes (25 July 2011—2011EQ, 2013EQ, and 2014EQ). The range of the vertical axes is adjusted to better visualize the majority of the data. Arrows indicate perturbations that occur between 2013EQ and 2014EQ. **B** Pore-pressure time series measured at depths below the seabed of 7.08 m (PZN-P6) and 6.28 m (PZS-P5). **C** Enlargement diagrams showing the pore-pressure evolutions during 2013EQ and 2014EQ. The horizontal scale uses the convention of day, month, and year.</figcaption></figure>

what was described previously³²⁻³⁵. Additionally, Fig. 4A highlights a 10-month period after the 2013EQ where the vertical displacement is almost locked with a displacement deficit of ±1 cm. Simultaneously and after the 2013EQ, pore pressure recorded by PZN-P6 synchronizes remarkably well with the TEKR GPS northern displacement indicating a strong link between the pore pressure from the MV (e.g., measured at PZN-P6) and the regional tectonic activity. The northern displacement shifts by around 0.35 cm to the north of the linear general trend that precedes 2013EQ (Fig. 2A).

Several authors have suggested, based on indirect evidence, a relationship between fluid activities and tectonic stresses⁵,¹⁷,²⁰,³⁶. Our data provide clear confirmation that regional strain (TEKR GPS data) and near-surface sediment pore pressure (PZN-P6) are interconnected (Fig. 4B). It is evident that a shallow localized process at the level of the piezometer is impossible to be detected

by onshore GPS data and vice versa. Therefore, the data suggest that shearing along the fault disturbs both the near-field sediment pore pressure and the far-field displacement.

The positive pore pressure recorded by PZN-P6 coincides with a displacement at TEKR toward the South (until 20/02/2014) while the pore-pressure decrease matches with a displacement to the North. In parallel, the locked vertical displacement at TEKR station during the wet period (high precipitations) indicates a net upward movement while the locked vertical displacement during the dry period implies a net downward vertical movement (Fig. 4A). The similarity between shallow pore water pressure and geodetic data suggests that pore-pressure perturbations have a physical and mechanical significance related to MMF at depth. We propose that this is due to the MV, which acts as a window to the MMF seismogenic zone linking stress/strain changes at depth to shallow pore-pressure variations.
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Pore-pressure versus geodetic data. **A** Vertical displacement from TEKR geodetic station (non-corrected for seasonal fluctuations—blue line) compared to Tekirdag daily mean temperatures (black squares) suggesting a seasonal ground surface displacement except for a 10-month period that begins in November 2013 (gray area). **B** Pore pressure measured at 7.08 mbsf during a 13-month period at PZN station (PZN-P6) synchronize after the 2013EQ with the TEKR GPS northern displacement (TEKR-North) corrected for seasonal fluctuations.</figcaption></figure>

**Negative pore pressure and dilatancy.** To illuminate the extent of the mechanical process controlling the near and far-field perturbations shown in Fig. 4, we hypothesize that the PZS-P5 negative pore pressure that developed only 1 week after the 2013EQ may be the result of an aseismic creeping activity along the MMF, which occurred in response to dilatancy due to shearing. Indeed, the PZS-P5 sensor is positioned within a silty-sandy layer (Supplementary Fig. 9) potentially dilatant during shearing. This shallow observed dilatancy is most probably the expression of a deeper and more generalized creeping process¹⁶, which affects the MV pore pressure and the regional displacement field. In contrast, the other five sensors at PZS positioned within clayey sediments did not show any pore-pressure perturbation during the monitoring period (Supplementary Fig. 9).

**Slow slip along the offshore fault disturbing the onshore geodetic data.** To test the shear-dilatancy hypothesis, we use a numerical model to describe the SSE. Using a purely elastic half-space dislocation model³⁷, we determine how displacement along the MMF affects surface deformation at TERK. To do this, the MMF was discretized into rectangles no larger than 2 km by 2 km. The deformation caused by different dislocations source types (i.e., pure right-lateral slip and a mix of right-lateral and dilatant slip) at different depths and locations were calculated (Methods and Supplementary Figs. 12–16). In order to reproduce the same trend as that observed at TEKR, the best fit is obtained with a slipping surface comprising of a shallow slipping patch (i.e., <8 km)

containing both right-lateral slip and a transient inflation component that propagates at 0.5 km/day from east to west (Fig. 5A and case C in Supplementary Table 1). These best fitting scenarios of a shallow slipping event are consistent with results from the literature obtained from the same region¹⁶,³⁸,³⁹. The model, however, predicts an east-west positive displacement not detectable on the TEKR data (Supplementary Fig. 2) suggesting that an isotropic and homogeneous elastic three-dimensional (3D) model with uniform slip may not be able to reproduce all the complexity of a natural system. The reality is more complex with a spatial variable slip; for example, less slip along the eastern section of MMF relative to TEKR compared to the western section would possibly reduce the simulated positive east-west motion at TEKR. An alternative explanation could be the amplification of north-south motion by normal faults to the north of MMF. These scenarios would require more complex models and data for validation.

By correlating the pore-pressure data from PZS-P5 with the longitude position of slipping (Fig. 5B), it is found that the 50-day recorded dilatancy occurred during a slipping event along the WH segment of the MMF (S1 in Fig. 5C). The shearing dilatancy and associated pore pressure increase the normal effective stress leading to a more stable slipping event⁴⁰. As a result, pore pressure slowly returns to equilibrium when the slow slip propagates to the west of WH. Proctor and co-authors⁴¹ have measured pore-pressure transients in hydraulically isolated laboratory faults during seismic and SSEs. Those authors showed the importance of pore fluid-rock interactions, which may dominate and control both stability and failure time of faults. Our field data confirm the
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>C</figcaption></figure>

<figure><img src="image_9.png" /></figure>

role of such pore fluid-sediment interactions in controlling or accompanying the process of fault slipping.

Discussion

Taken all those observed and modeling results together, it is possible to draw a timeline of events that occurred after the

2013EQ indicating the way the aseismic creep affects the MV activities as well as the 3D displacement field surrounding the MMF (Fig. 5C). Shallow slipping along the MMF (within the syn- and pre-kinematic layers42) propagating from east to west is probably at the origin of the far-field observed geodetic data from TEKR station: up and south displacements before 20/02/2014 and down and north displacements after this date (Fig. 5C). The
**Fig. 5** Time line scenario. A Calculated north displacement at TEKR versus longitude position of slipping compared to observed north displacement versus time. Yellow curves correspond to the displacement towards the south while the blue curves correspond to the displacement toward the north. Pore-pressure data from PZN-P6 versus time are also shown with an increase in the pore pressure during the first stage corresponding to the increase of the normal stress at the outer limit of the damaged zone (yellow dots). Earthquake (2013EQ) timing is indicated by a triangle (top horizontal axis) and its position is marked by a black rectangle (base horizontal axis). The TEKR geodetic station is indicated by a black rectangle (base horizontal axis). B Calculated up displacement at TEKR (yellow corresponds to up and blue to down) versus longitude position of slipping compared to pore pressure from piezometer PZS-P5 with a decrease of the pore pressure during the slip dilatancy affecting the piezometer area (yellow dots). A 50-day dilatancy period corresponds to a creep dilatancy along the segment S1. C 3D diagram illustrating the consequence of aseismic creep of the strike-slip fault on the mud volcano activities and the 3D displacement field surrounding the Main Marmara Fault (MMF). A black star indicates the 2013EQ earthquake epicenter and a blue dot the TEKR geodetic station. Locations of piezometers PZN and PZS are also indicated. The aseismic slipping propagates along the MMF at 0.5 km/day from east to west lasting >10 months.

shear-dilatancy caused by the slow slipping along the WH fault segment recorded by PZN-P5 between 2/12/2013 and 20/01/2014 corresponds to the reactivation of the MV by increasing the normal stress at the outer limit of the damaged zone. The MV pressurized by the increase of the normal stress is at the origin of the high pore pressure recorded by PZN-P6 (red dots in Fig. 5A). The subsequent decay of the pore pressure recorded by PZN-P6 (blue dots in Fig. 5A) reaching negative values is most likely the result of the decrease of the normal stress at the MV boundary causing the swelling of the MV and requiring the replacement of the dissipated pore-fluid volume.

This aseismic slipping event propagates along the MMF at 0.5 km/day from east to west lasting more than 10 months. However, based solely on the available data, it is not possible to ascertain when or where this 10-month SSE has nucleated. It may have started with the 2013EQ or alternatively the 2013 earthquake may have been triggered by the SSE. The slipping event we observe occurs mainly along the WH segment, a region having the largest proportion of mainshocks with associated foreshocks and aftershocks<sup>43</sup>. According to Martinez-Garzon and co-authors<sup>43</sup>, the WH segment is close to failure increasing therefore the susceptibility of the segment to seismic triggering. Our results confirm also what was observed by Caniven and co-authors<sup>44</sup>, indicating based on numerical calculations that fault dilation is a "key factor controlling fault slip stability."

Comparable features have been observed for SSE on subduction zones<sup>45</sup>. Shorter duration (3–20 days) SSEs have been detected along the San Andreas fault, however only at depth (>16 km)<sup>46</sup>. It has been proposed that SSEs have been linked with the migration of pore-pressure waves in subduction zones<sup>47</sup> and strike-slip faults<sup>48</sup> caused by fault valving. The transient nature of the observed inflation at PZN and PZS as well as the geodetic displacement obtained from a single GPS station (TEKR) provide direct evidence supporting this proposition.

Beyond the present case concerning the MMF, our data provide valuable information on the kinematics of transform faults. The deployment of surface piezometers measuring pore pressure all along a transform fault can help define the time scale regulating the coupling between faulting, creeping, and dilatancy processes. MVs are common features affecting the seafloor, particularly at subduction zones, where the relationship between MV plumbing systems and slow-slip mechanisms is an emerging field in earthquake research. MV instrumentation with surface piezometers may provide valuable information on understanding the shear-strain coupling at depth, within the seismogenic zone. The pressure drop recorded 4 days before the 2013EQ suggests that such an experiment could be a source of improving scientific knowledge in earthquake nucleation processes. We also show that new insights can be gained from the detailed observation of onshore geodetic data since small disturbances seem to have major meaning in terms of transient deformations along active faults.

## Methods

**Seabed amplitude and sub-seabed seismic features.** Site surveys including high-resolution 3D seismic and AUV bathymetric data were used<sup>15</sup> to draw the seabed amplitude and sub-seabed seismic features characterizing the study area (Supplementary Fig. 1).

**Geodetic data.** The 7-year period GPS horizontal displacement (east with respect to fixed EU plate) time series<sup>6</sup> shows that the Anatolian plate evolves smoothly to the west relatively independently of 2013EQ (Supplementary Fig. 2). Horizontal displacements (east) recorded by the two stations from the Eurasian plate (ISTN and TEKR) remain almost unchanged during the monitoring period. The 2014EQ unequally affects data from BALK, YENC, and CANA stations.

**Tools: piezometer and coring.** Free-fall piezometer equipped with a sediment-piercing lance of 60 mm diameter<sup>49</sup> is used to perform in situ pore pressure and temperature measurements. The piezometer is ballasted with lead weights (up to 1000 kg) and equipped with differential pore pressure and temperature sensors. Two piezometers were deployed within the shallow sub-surface (<7.08 mbsf) at sites PZN and PZS. The piezometer pore pressure and temperature sensors have an accuracy of ±0.5 kPa and 0.05 °C, respectively (data in Supplementary Figs. 3–8).

Core MNT-KS26 was recovered at 300 m to the west of PZS (Fig. 1). The total length of the core is around 8 m. Density and P-wave velocity data were obtained on undisturbed sediment sections using the GEOTEK core logging devices (MSCL). Supplementary Fig. 9a, b shows the density and P-wave velocity versus depth. Intervals characterized by high density and/or high P-wave velocity anomalies are indicated as gray areas and could correspond to coarser materials. The effective lithostatic stress profile shown in Supplementary Fig. 9c was calculated from the density profile.

**Pore-pressure and temperature data.** Pore-pressure and temperature data from PZN recorded for more than 1 year are shown in Supplementary Fig. 3. The pore pressure generated by the piezometer penetration decays smoothly during several weeks and stabilizes for the shallowest four sensors (Supplementary Fig. 3a). Pore-pressure data from the two deepest two sensors fluctuate throughout the monitoring period (Supplementary Fig. 3a). Negative pore pressures as low as −13 kPa were recorded by the deepest sensor (Supplementary Fig. 3a). The temperature decays rapidly to reach in situ equilibrium temperature (Supplementary Fig. 3b). The temperatures at the deepest three sensors intersect regularly and are almost constant throughout the monitoring period.

Supplementary Fig. 4 shows an enlargement of pore-pressure generation during 2013EQ (Supplementary Fig. 4a) and 2014EQ (Supplementary Fig. 4b). During 2013EQ, pore-pressure data from the 6 sensors were affected for a few hours and fluctuates between −4 and +4 kPa. Despite the large distance of PZN from the 2014EQ epicenter (Fig. 1), the earthquake significantly affects the pore pressure from the 6 sensors which fluctuate between −25 and +8 kPa (Supplementary Fig. 4b).

Pore-pressure and temperature data recorded for >1 year at PZS are shown in Supplementary Fig. 5. The pore pressure and temperature generated by the piezometer penetration decays rapidly after installation (Supplementary Fig. 5a). Pore-pressure data from P5 fluctuate between −8 and 0 kPa throughout the monitoring period (Supplementary Fig. 5a).

Supplementary Fig. 6 shows an enlargement of pore-pressure generation during 2013EQ (Supplementary Fig. 6a) and 2014EQ (Supplementary Fig. 6b). Pore pressure data from the 6 sensors were affected by both events and fluctuates between −1 and +1 kPa during 2013EQ (Supplementary Fig. 6a) and between −7 and +8 kPa during 2014EQ (Supplementary Fig. 6b).

Thermal and hydraulic gradients plotted in Supplementary Fig. 7 correspond to ten different periods of monitoring. Thermal gradients at PZN (Supplementary Fig. 7a) suggest that the temperature field at the level of the upper four sensors is in a permanent regime and primarily diffusion-controlled (quasi-constant gradient) while it is advection-controlled along the deepest two sensors (quasi-constant temperature). Hydraulic gradients shown in Supplementary Fig. 7a indicate that
the pore pressure perturbation only concerns the two deepest sensors. Thermal data from PZS (Supplementary Fig. 7b) indicate a transient temperature regime without the possibility to conclude about the thermal process (advection or diffusion) controlling the temperature field. The non-linearity of the temperature profile could be the result of a transient diffusion of advection processes. The hydraulic gradients at PZS show that only the sensor P5 at 6.28 mbsf is concerned by the pore pressure perturbations.

The measured initial pore water pressure generated by the rod insertion ($\Delta u_i$) depends on the shear strength and the elastic properties of the medium$^{50}$. Data in Supplementary Fig. 8a show a general linear increase of $\Delta u_i$ with depth except for PZS-P5 and PZN-P6 deviating slightly from this trend and indicating stiffer materials. The hydraulic diffusivity $C_h$ (or the horizontal coefficient of consolidation) of the medium normalized by the square root of the rigidity index ($\sqrt{I_r}$) was derived by using the cavity expansion theory$^{49}$ and by back-analyzing the pore-pressure decay curves following the piezometer penetration. Comparisons between the measured and the calculated pore pressure decay at PZN and PZS are shown, respectively, in Supplementary Fig. 8c, d.

**Transient diffusion–advection process.** To understand how pore-pressure variations at PZN-P6 and fluid advection may affect the pore-pressure field in the above sedimentary layers, a one-dimensional (1D) transient diffusion–advection equation (Eq. 1) is solved numerically using a finite difference method.

$$ \frac{\partial u}{\partial t} = \frac{\partial}{\partial z} \left( C \frac{\partial u}{\partial z} \right) + v_z \frac{\partial u}{\partial z} \quad (1) $$

In Eq. (1), C signifies the hydraulic diffusivity of the sediment, $v_z$ the vertical velocity of the fluid, $u$ the pressure field, $t$ the simulation time, and $z$ the depth below the seabed. To solve numerically the 1D diffusion–advection equation, a centered explicit finite difference discretization scheme is used. In this calculation, C was taken equal to $1.5 \times 10^{-8} \text{ m}^2/\text{s}$ equivalent to the $C_h$ values at PZN-P6 (Supplementary Fig. 8b) and for a rigidity index $I_r$ of 40. Three values of $v_z$ are considered showing that pure diffusion ($v_z = 0 \text{ m/s}$) is needed to reproduce the general trend of PZN-P5 data while $v_z$ value of $3 \times 10^{-8} \text{ m/s}$ is required to simulate the observed peaks at PZN-P5. For the considered C and $v_z$ values, the PZN-P6 imposed pressure affected none of the levels above PZN-P5. Those results confirm that the hydraulic process controlling the pore-pressure field above PZN-P6 is an alternation between pure diffusion and diffusion–advection and this depends on the pressure level and hydro-fracturing at PZN-P6.

**Pore pressure and hydro-fracturing.** Hydro-fracturing in sediments takes place when the excess pore pressure exceeds the least principle stress plus the tensile strength of the medium. Comparing recorded pore pressure at PZN-P5 (Supplementary Fig. 11a) and PZN-P6 (Supplementary Fig. 11b) with effective lithostatic stresses ($\sigma'_v$) at corresponding depth indicates that PZN-P6 pore-pressure peaks are reaching 70–100% of $\sigma'_v$ before dropping sharply. These high pore-pressure values suggest that hydro-fracturing occurred at PZN-P6 level during the first 5 months of monitoring. In contrast, the pore pressure at PZN-P5 is too low to initiate hydro-fracturing.

**3D displacement field in the half-space linear elastic medium due to shear and tensile along the MMF.** The aim of this modeling was to ascertain whether the anomalous deformation observed at the geodetic station TEKR could be explained by a slipping event on the MMF and if so, provide some constraint on it. To do this, a 3D elastic half-space dislocation model$^{37}$ with a Poisson ratio of 0.25 was used to model the surface displacement. This model is based on elastic dislocation theory whereby the displacement in the i direction at position x due to a dislocation d on a planar fault subdivided into N individual cells $\Gamma(i)$:

$$ \Delta u_i(x) = \sum_{i=1}^{N} \left( \int_{\Gamma(i)} T_{ik}^{ss}(x, \xi_i) n_k d_{ss}(\xi_i) dS + \int_{\Gamma(i)} T_{ik}^{t}(x, \xi_i) n_k d_i(\xi_i) dS \right) \quad (2) $$

where $n_k$ is the normal vector to the fault plane and $d_{ss}$ and $d_i$ are shear and tensile slip in cell i, respectively. $T_{ik}^{ss}(x, \xi_i)$ and $T_{ik}^{t}(x, \xi_i)$ map the shear and tensile slip at position $\xi$ to the displacement in the i direction at position x$^{51}$ using the method of images with analytical description of deformation in an elastic medium, these kernels have explicit formulas for both shear and tensile rectangular sources$^{52}$. The principle of superposition allows different source mechanisms and slipping surfaces to be calculated individually and summed together to produce the total displacement. As a result, the digitized fault trace$^{11}$ of the MMF is subdivided into segments no larger than 2 km. Assuming the MMF has a dip of 90°, a discretized fault plane was constructed with individual slipping cells no larger than 2 km × 2 km producing $158 \times 7$ non-overlapping cells. Two faulting mechanisms are considered (a) permanent pure strike slip (i.e., $d_t = 0$) and (b) a permanent strike slip with a transitory tensile component in the active cell (i.e., $d_i$ is a box-car function in time while $d_{ss}$ is a Heaviside function). The choice of pure-right lateral slip is in keeping with the regional geodetic motion$^{31}$ as well as the historical seismicity$^{14}$. The inflation contributions from the different elements are not summated as we assume, based on observations at PZN and PZS, that inflation is a transient phenomenon only active during slipping.

To investigate the propagation of a slipping front across the fault, the contribution from the cells is summed sequentially based on geographical location (i.e., from east to west and from west to east) allowing for the assessment that different sections of the fault make the displacement at TEKR. Therefore to ascertain the displacement in the direction i at TEKR due to slipping event propagating from east to west where the active slipping cell is cell n, Eq. 2 can be rewritten as:

$$ u_i^{EW}(n) = r \sum_{i=1}^{n \le N} \Delta u_i^s + (r-1) \Delta u_i^t \quad (3) $$

where $\Delta u_i^s$ is the displacement due to 1 m strike slip of cell i (i.e., the first integral in Eq. (2)), and $\Delta u_i^t$ is the contribution from the 1 m of tensile displacement of cell i with the cells ordered sequentially from east to west. r represents the ratio of strike to tensile slip. Using the same formalism, west to east deformation is defined as:

$$ u_i^{WE}(n) = r \sum_{i=1}^{n \le N} \Delta u_{N+1-i}^s + (r-1) \Delta u_{N+1-i}^t \quad (4) $$

Sensitivity to depth was tested by propagating the slipping event across the fault at a constant depth for each simulation with a width of 1 cell (i.e., 2 km). The procedure for west-to-east simulation is depicted in Supplementary Fig. 12. In summary, four cases are considered, Case A; a purely right-lateral rupture propagates from east to west; Case B: a purely right-lateral rupture propagates from west to east; Case C: a right-lateral rupture with a tensile component propagates from east to west; Case D: the right-lateral rupture with a tensile component propagates from west to east. Within each of these cases, the slip surface has been placed at a range of different depths.

With only one observation point, we have focused our analysis on comparing the simulated general trends of the vertical and north-south displacement at the location of TEKR with the observations. The key features we are looking to observe are an initial displacement to the north followed by a motion to the south, and on the vertical component, an upward displacement followed by subsidence. Supplementary Table 1 summarizes all the simulations and whether they conform to the aforementioned trends while Supplementary Figs. 13–16 depict all simulated ground motion at TEKR. As shown in Supplementary Table 1, the only simulations that match these trends is the case where the rupture propagates from east to west and rupture occurs at a shallow depth (i.e., 0–6 km in the pure strike slip case, and 0–8 km in the mixed strike slip–tensile rupture).

In order to produce the observed southern displacement of 0.35 cm at TEKR, this would require a right-lateral and dilatant displacement of 0.6 cm; and 0.2 cm in the case of a pure strike-slip movement, which are equivalent to M5.3 and M5.1 earthquakes, respectively (using a shear modulus of 30 GPa and considering them as one large continuous event). However, changing the element size (e.g., increasing/decreasing cell size, including multiple cells from different depths with different slips at the same time) or the inclusion of normal faults in the vicinity of MMF may change these values. Consequently, no conclusion is drawn on the size of the active slipping zone.

The distance between the barycenter of the sources is converted to time by assuming a constant rupture velocity for the slipping event, the value that provided the best fit to the observations was a velocity of 0.5 km/day.

## Data availability

GPS data are freely available at Nevada Geodetic Laboratory (http://geodesy.unr.edu) and piezometer data can be downloaded from https://doi.org/10.17882/79781. Geological (AUV bathymetry, seismic lines) and sedimentological data are available in the main text or the Supplementary Materials.

## Code availability

Computer code used to generate results that are reported in Supplementary Figs. 8, 10, and 13–16 are available at https://github.com/nsultan-2021/advection-diffusion. Surface displacement calculations can be found here: https://github.com/s-murfy/StrikeSlipDef

Received: 16 June 2021; Accepted: 16 March 2022;

Published online: 08 April 2022

## References

1. Dragert, H., Wang, K. & James, T. S. A silent slip event on the deeper Cascadia subduction interface. *Science* **292**, 1525–1528 (2001).

2. Kanamori, H. & Brodsky, E. E. The physics of earthquakes. *Rep. Prog. Phys.* **67**, 1429 (2004).

3. Segall, P., Rubin, A. M., Bradley, A. M. & Rice, J. R. Dilatant strengthening as a mechanism for slow slip events. *J. Geophys. Res. Solid Earth* https://doi.org/10.1029/2010JB007449 (2010).

4. Ikari, M. J., Marone, C., Saffer, D. M. & Kopf, A. J. Slip weakening as a mechanism for slow earthquakes. *Nat. Geosci.* **6**, 468–472 (2013).
5. Jónsson, S., Segall, P., Pedersen, R. & Björnsson, G. Post-earthquake ground movements correlated to pore-pressure transients. *Nature* **424**, 179–183 (2003).

6. Blewitt, G., Hammond, W. C. & Kreemer, C. Harnessing the GPS data explosion for interdisciplinary science. *Eos* **99**, 1–2 (2018).

7. Obara, K. Nonvolcanic deep tremor associated with subduction in southwest Japan. *Science* **296**, 1679–1681 (2002).

8. Gosselin, J. M. et al. Seismic evidence for megathrust fault-valve behavior during episodic tremor and slip. *Sci. Adv.* **6**, eaay5174 (2020).

9. Warren-Smith, E. et al. Episodic stress and fluid pressure cycling in subducting oceanic crust during slow slip. *Nat. Geosci.* **12**, 475–481 (2019).

10. Le Pichon, X. et al. The active main Marmara fault. *Earth Planet. Sci. Lett.* **192**, 595–616 (2001).

11. Şengör, A. M. C. et al. The geometry of the North Anatolian transform fault in the Sea of Marmara and its temporal evolution: implications for the development of intracontinental transform faults. *Can. J. Earth Sci.* **51**, 222–242 (2014).

12. Rangin, C., Le Pichon, X., Demirbag, E. & Imren, C. Strain localization in the Sea of Marmara: propagation of the North Anatolian Fault in a now inactive pull-apart. *Tectonics* https://doi.org/10.1029/2002TC001437 (2004).

13. Géli, L., Henry, P. & Çağatay, M. N. A review of 20 years (1999–2019) of Turkish–French collaboration in marine geoscience research in the Sea of Marmara. *Med. Geosci.* *Rev.* **3**, 3–27 (2021).

14. Schmittbuhl, J., Karabulut, H., Lengliné, O. & Bouchon, M. Seismicity distribution and locking depth along the Main Marmara Fault, Turkey. *Geochem. Geophysi. Geosyst.* **17**, 954–965 (2016).

15. Uchida, N., Kalafat, D., Pinar, A. & Yamamoto, Y. Repeating earthquakes and interplate coupling along the western part of the North Anatolian Fault. *Tectonophysics* **769**, 228185 (2019).

16. Yamamoto, R. et al. Seafloor geodesy revealed partial creep of the North Anatolian Fault submerged in the Sea of Marmara. *Geophys. Res. Lett.* **46**, 1268–1275 (2019).

17. Dupre, S. et al. Tectonic and sedimentary controls on widespread gas emissions in the Sea of Marmara: results from systematic, shipborne multibeam echo sounder water column imaging. *J. Geophys. Res. Solid Earth* **120**, 2891–2912 (2015).

18. Thomas, Y. et al. Contribution of high-resolution 3D seismic near-seafloor imaging to reservoir-scale studies: application to the active North Anatolian Fault, Sea of Marmara. *Surf. Geophys.* **10**, 291–301 (2012).

19. Ruffine, L. et al. Multiple gas reservoirs are responsible for the gas emissions along the Marmara fault network. *Deep Sea Res. Part II Top. Stud. Oceanogr.* **153**, 48–60 (2018).

20. Henry, P. et al. A statistical approach to relationships between fluid emissions and faults: the Sea of Marmara case. *Deep Sea Res. Part II Top. Stud. Oceanogr.* **153**, 131–143 (2018).

21. Geli, L. et al. Gas and seismicity within the Istanbul seismic gap. *Sci. Rep.* https://doi.org/10.1038/s41598-018-23536-7 (2018).

22. Sultan, N., Plaza-Faverola, A., Vadakkekuliyambatta, S., Buenz, S. & Knies, J. Impact of tides and sea-level on deep-sea Arctic methane emissions. *Nat. Commun.* **11**, 1–10 (2020).

23. Geli, L., Henry, P., Çağatay, N., Apprioual, R. & Mickael, R. Sea Open Scientific Data Publication - SEANOE (2021). https://doi.org/10.17882/79781.

24. Wollin, C., Bohnhoff, M., Martínez-Garzón, P., Küperkoch, L. & Raub, C. A unified earthquake catalogue for the Sea of Marmara Region, Turkey, based on automatized phase picking and travel-time inversion: seismotectonic implications. *Tectonophysics* **747**, 416–444 (2018).

25. Saltogianni, V., Gianniou, M., Taymaz, T., Yolsal-Çevikbilen, S. & Stiros, S. Fault slip source models for the 2014 Mw 6.9 Samothraki-Gökçeada earthquake (North Aegean trough) combining geodetic and seismological observations. *J. Geophys. Res. Solid Earth* **120**, 8610–8622 (2015).

26. Mekik, C., Yildirim, O. & Bakici, S. The Turkish real time kinematic GPS network (TUSAGA-Aktif) infrastructure. *Sci. Res. Essays* **6**, 3986–3999 (2011).

27. Altamimi, Z., Rebischung, P., Métivier, L. & Collilieux, X. ITRF2014: a new release of the International Terrestrial Reference Frame modeling nonlinear station motions. *J. Geophys. Res. Solid Earth* **121**, 6109–6131 (2016).

28. Doğan, U. Accuracy analysis of relative positions of permanent GPS stations in the Marmara region, Turkey. *Surv. Rev.* **39**, 156–165 (2007).

29. Cleveland, R. B., Cleveland, W. S., McRae, J. E. & Terpenning, I. STL: a seasonal-trend decomposition. *J. Off. Stat.* **6**, 3–73 (1990).

30. R Core Team. *R: A Language and Environment for Statistical Computing* (R Foundation for Statistical Computing, 2013).

31. England, P., Houseman, G. & Nocquet, J. M. Constraints from GPS measurements on the dynamics of deformation in Anatolia and the Aegean. *J. Geophys. Res. Solid Earth* **121**, 8888–8916 (2016).

32. Amos, C. B. et al. Uplift and seismicity driven by groundwater depletion in central California. *Nature* **509**, 483–486 (2014).

33. Blewitt, G., Lavallée, D., Clarke, P. & Nurutdinov, K. A new global mode of Earth deformation: Seasonal cycle detected. *Science* **294**, 2342–2345 (2001).

34. Johnson, C. W., Fu, Y. & Bürgmann, R. Seasonal water storage, stress modulation, and California seismicity. *Science* **356**, 1161–1164 (2017).

35. Larochelle, S., Gualandi, A., Chanard, K. & Avouac, J. P. Identification and extraction of seasonal geodetic signals due to surface load variations. *J. Geophys. Res. Solid Earth* **123**, 11–031 (2018).

36. Plaza-Faverola, A. & Keiding, M. Correlation between tectonic stress regimes and methane seepage on the western Svalbard margin. *Solid Earth* **10**, 79–94 (2019).

37. Okada, Y. Surface deformation due to shear and tensile faults in a half-space. *Bull. Seismol. Soc. Am.* **75**, 1135–1154 (1985).

38. Rousset, B. et al. An aseismic slip transient on the North Anatolian Fault. *Geophys. Res. Lett.* **43**, 3254–3262 (2016).

39. Aslan, G. et al. Shallow creep along the 1999 Izmit Earthquake rupture (Turkey) from GPS and high temporal resolution interferometric synthetic aperture radar data (2011–2017). *J. Geophys. Res. Solid Earth* **124**, 2218–2236 (2019).

40. Segall, P. & Rice, J. R. Dilatancy, compaction, and slip instability of a fluid-infiltrated fault. *J. Geophys. Res. Solid Earth* **100**, 22155–22171 (1995).

41. Proctor, B., Lockner, D. A., Kilgore, B. D., Mitchell, T. M. & Beeler, N. M. Direct evidence for fluid pressure, dilatancy, and compaction affecting slip in isolated faults. *Geophys. Res. Lett.* **47**, e2019GL086767 (2020).

42. Kende, J., Henry, P., Bayrakci, G., Özeren, M. S. & Grall, C. Moho depth and crustal thinning in the Marmara Sea region from gravity data inversion. *J. Geophys. Res. Solid Earth* **122**, 1381–1401 (2017).

43. Martínez-Garzón, P., Ben-Zion, Y., Zaliapin, I. & Bohnhoff, M. Seismic clustering in the Sea of Marmara: implications for monitoring earthquake processes. *Tectonophysics* **768**, 228176 (2019).

44. Caniven, Y., Morgan, J. K. & Blank, D. G. The role of along-fault dilatancy in fault slip behavior. *J. Geophys. Res. Solid Earth* **126**, e2021JB022310 (2021).

45. Ragu diet, M. et al. Triggering of the 2014 M w 7.3 Papanoa earthquake by a slow slip event in Guerrero, Mexico. *Nat. Geosci.* **9**, 829–833 (2016).

46. Rousset, B., Bürgmann, R. & Campillo, M. Slow slip events in the roots of the San Andreas fault. *Sci. Adv.* **5**, eaav3274 (2019).

47. Cruz-Atienza, V. M., Villafuerte, C. & Bhat, H. S. Rapid tremor migration and pore-pressure waves in subduction zones. *Nat. Commun.* **9**, 1–13 (2018).

48. Zhu, Y., Wang, K. & He, J. Effects of earthquake recurrence on localization of interseismic deformation around locked strike-slip faults. *J. Geophys. Res. Solid Earth* **125**, e2020JB019817 (2020).

49. Sultan, N. & Lafuerza, S. In situ equilibrium pore-water pressures derived from partial piezoprobe dissipation tests in marine sediments. *Can. Geotech.* **J.** **50**, 1294–1305 (2013).

50. Burns, S. E. & Mayne, P. W. Pore pressure dissipation behavior surrounding driven piles and cone penetrometers. *Transportation Res. Rec.* **1675**, 17–23 (1999).

51. Gimbutas, Z., Greengard, L., Barall, M. & Tullis, T. E. On the calculation of displacement, stress, and strain induced by triangular dislocations. *Bull. Seismol. Soc. Am.* **102**, 2776–2780 (2012).

52. Okada, Y. Internal deformation due to shear and tensile faults in a half-space. *Bull. Seismol. Soc. Am.* **82**, 1018–1040 (1992).

53. Uieda, L. et al. PyGMT: a Python interface for the generic mapping tools (PyGMT, 2021).

Acknowledgements

The authors thank Nurcan Meral Özel (KOERI), P. I. of MARSITE project. Luca Gasperini (CNR), led the deployment cruise (MARMARA-2013) of R/V Urania. Livio Ruffine, Namik Çagatay, and Pierre Henry co-led the recovery cruise (MARSITECRUISE) of R/V Pourquoi pas? Bruno Marsset, Yannick Thomas, and Stephan Ker are thanked for providing the 3D seismic data. The marine operations for instrument deployment, recovery, and analyses were supported by the 7th EU-funded Framework Programme (FP7), within EMSO and MARSITE (EU Grant agreement ID: 308417) Projects and ANR MAREGAMI project (Grant No. ANR-16-CE03-0010). Supplementary Fig. 12 is generated with PyGMT<sup>53</sup>.

Author contributions

N.S. led the piezometer and geodetic data analyses and wrote the manuscript. S.M. carried out the 3D numerical calculations to derive the 3D displacement field at TEKR and contributed to data analyses and interpretations. V.R. carried out the geological analyses and contributed to data interpretations. L.G. conceived and designed the off-shore experiment and contributed to data analyses and interpretations. S.M., V.R., and L.G. helped write the manuscript.

Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-022-29558-0.

**Correspondence** and requests for materials should be addressed to Nabil Sultan.

**Peer review information** *Nature Communications* thanks Sinan Ozeren and the other anonymous reviewer(s) for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022