Elastic Crust Research Articles

Characterizing the physical properties of gelatin, a classic analog for the brittle elastic crust, insight from numerical modeling

Precise characterization of the mechanical properties of gelatin, a classic analog of the elastic crust, is necessary for scaling the mechanical models of the Earth's crust behavior in laboratory experiments. Here we reassess how to accurately calculate the Young modulus (E) of gelatin contained in experimental tanks. By means of dedicated analog experiments and finite element simulations, we estimate the bias introduced by using equations appropriate for a half-space to interpret the subsidence due to a cylindrical surface load applied on the gelatin. In the case of a standard experimental setup with gelatin adhering to the tank wall, we find E is overestimated by at least 5% for a box with lateral size smaller than 20 times the cylinder diameter. In addition, we deduce a correction factor to be applied when using an analytical formula. We confirm that measuring the shear velocity leads to accurate estimates for the rigidity of gelatin. We also propose a new method for in situ Young's modulus estimation, relying on the length of air-filled propagating crack. Indeed, for a given injected volume, this length depends only on the density contrast between air and gelatin and on the Young's modulus of the gelatin. The fracture toughness of the gelatin is estimated independently. Direct comparison between fracture toughness and Young's modulus shows that for a given Young's modulus, salted gelatin has a higher fracture toughness than unsalted gelatin.

Dynamical tides in neutron stars: the impact of the crust

ABSTRACT We consider the dynamical tidal response of a neutron star in an inspiralling binary, focussing on the impact of the star’s elastic crust. Within the context of Newtonian gravity, we add the elastic aspects to the theoretical formulation of the problem and quantify the dynamical excitation of different classes of oscillation modes. The results demonstrate the expectation that the fundamental mode dominates the tidal response and show how the usual tidal deformability (and the Love number) emerge in the static limit. In addition, we consider to what extent the different modes may be excited to a level where the breaking strain of the crust would be exceeded (locally). The results show that the fundamental mode may fracture the crust during the late stages of inspiral. This is also the case for the first gravity mode, which reaches the breaking threshold in strongly stratified stars. In our models with a fluid ocean, interface modes associated with the crust-ocean transition may also induce crust fracture. If this happens, it does so earlier in the inspiral, at a lower orbital frequency.

A Superfluid Perspective on Neutron Star Dynamics

As mature neutron stars are cold (on the relevant temperature scale), one has to carefully consider the state of matter in their interior. The outer kilometre or so is expected to freeze to form an elastic crust of increasingly neutron-rich nuclei, coexisting with a superfluid neutron component, while the star’s fluid core contains a mixed superfluid/superconductor. The dynamics of the star depend heavily on the parameters associated with the different phases. The presence of superfluidity brings new degrees of freedom—in essence we are dealing with a complex multi-fluid system—and additional features: bulk rotation is supported by a dense array of quantised vortices, which introduce dissipation via mutual friction, and the motion of the superfluid is affected by the so-called entrainment effect. This brief survey provides an introduction to—along with a commentary on our current understanding of—these dynamical aspects, paying particular attention to the role of entrainment, and outlines the impact of superfluidity on neutron-star seismology.

Present‐Day and Long‐Term Uplift Across the Western Transverse Ranges of Southern California

AbstractIt has been known for decades that the present‐day shortening rates across the Western Transverse Ranges (WTR) in Southern California are rapid, reaching 10–15 mm/year near the heavily populated Los Angeles area. However, only recently have geodetic measurements of vertical motion in the WTR been sufficiently dense to resolve a tectonic vertical signal. In this study, we show that much of the geodetically derived vertical velocity field in the WTR can be attributed to the interseismic signal of strain accumulation on reverse faults. We invert geodetic and geologic data for slip rate and interseismic coupling on faults using a kinematic model consisting of faults embedded in an elastic crust over an inviscid mantle. This method allows us to infer the permanent, long‐term component of vertical motions from recoverable, short‐term motions. We infer that much of the geodetically observed 3–4 mm/year of differential vertical motion across the WTR, involving subsidence along the Santa Barbara coastline and uplift of the Santa Ynez Range, can be attributed to recoverable elastic deformation associated with interseismic locking on faults dipping under the WTR. The sum of dip‐slip rates across the WTR decreases from 10.5–14.6 mm/year on the east side near Ventura, California to 5–6.2 mm/year across the western side of the Santa Barbara Channel. The total moment accumulation rate in both the Santa Barbara Channel and the combined San Fernando Valley‐LA Basin regions is equivalent to about two Mw = 7 earthquakes every 100 years.

Recent Movements of the Earth’s Surface and Stress Propagation in the Elastic Upper Crust

The Earth’s crust is modeled by a thin elastic plate, whereas the underlying lithosphere by a half-space with viscous fluid rheology. For this system, with the use of the Fourier transform with respect to horizontal coordinates and the Laplace transform with respect to time, the solutions of continuum mechanics equations are obtained in the form of the diffusion type waves which propagate with strong attenuation from the initial perturbation domain across the Earth’s surface and cause its displacements. Analytical formula for these waves is obtained for the case of a point initial perturbation, which explicitly links the Earth’s surface displacements and stresses in the elastic crust with the horizontal coordinates and time. The diffusion waves (inertialess Rayleigh and Love waves) can be considered as a mechanism of recent vertical and horizontal movements of the Earth’s crust. The inertialess Rayleigh waves restore the isostatic equilibrium of the Earth’s crust which is disturbed by the initial vertical displacement of the Earth’s surface. In the case when an inertialess Love wave propagates along a locked fault from a segment where the fault has ruptured which was accompanied by an earthquake, this wave lowers the normal stress applied to the fault’s side and, hence, reduces the friction force thus facilitating rupture on a remote part of the locked fault and causing a new earthquake.

Are normal fault earthquakes due to elastic rebound or gravitational collapse?

We discuss two competing models for explaining the ground deformation associated with normal faulting earthquake in the brittle elastic upper crust. The classic elastic rebound theory is usually applied for all tectonic settings. In normal fault earthquakes, this model would predict a horizontal stretching eventually responsible for the elastic rebound at the earthquake. However, volumes mostly subside vertically during an extensional earthquake and the collapsed ground in the hanging wall is about one order of magnitude larger than the uplifted volumes of the surrounding hanging wall and footwall. The elastic rebound model would explain this asymmetry with a high horizontal elastic compressibility of the hanging wall and footwall absorbing the coseismic push. We rather suggest that the force activating normal fault earthquakes is mostly dictated by the sliding of the hanging wall, owing gravitational potential. The much larger coseismic subsidence with respect to the uplift can be explained by the closure at depth of a diffuse network of microfractures developed during the interseismic period. Since the horizontal stretching does not exist below ~1 km of depth, with the minimum horizontal stress tensor becoming positive below that depth, the development of a normal fault can be activated only by the vertical maximum stress tensor, i.e., the lithostatic load. The common fluids expulsion at the coseismic stage requires diffuse secondary permeability in the upper crust, in agreement with the presence of a diffuse network of microfractures.

Tidal deformations of neutron stars with elastic crusts

With the first detections of binary neutron star mergers by gravitational-wave detectors, it proves timely to consider how the internal structure of neutron stars affects the way in which they can be asymmetrically deformed. Such deformations may leave measurable imprints on gravitational-wave signals and can be sourced through tidal interactions or the formation of mountains. We detail the formalism that describes fully-relativistic neutron star models with elastic crusts undergoing static perturbations. This formalism primes the problem for studies into a variety of mechanisms that can deform a neutron star. We present results for a barotropic equation of state and a realistic model for the elastic crust, which enables us to compute relevant quantities such as the tidal deformability parameter. We find that the inclusion of an elastic crust provides a very small correction to the tidal deformability. The results allow us to demonstrate when and where the crust starts to fail during a binary inspiral and we find that the majority of the crust will remain intact up until merger.

Impact of power-law rheology on the viscoelastic relaxation pattern and afterslip distribution following the 2010 Mw 8.8 Maule earthquake

After large earthquakes at subduction zones, the plate interface continues moving due to mostly frictional afterslip or simply afterslip processes. Below approximately 60 km depth, the seismic moment release at the plate interface is quite small indicating that the shear strength is low and stable sliding is the prevailing process. This agrees with the lack of significant interseismic locking at deeper segments (>60 km) resulting from the inversion of geodetic data and thus low afterslip can be expected. However, inversion models that employ linear viscoelastic mantle rheology and an elastic crust result in significant afterslip at depths >60 km. In this paper, we present a combination of a 3D forward geomechanical model with power-law rheology that simulates postseismic relaxation with dislocation creep processes in the crust and upper mantle and an afterslip inversion. We estimate the cumulative viscoelastic relaxation and the afterslip distribution for the first six years following the 2010 Mw 8.8 Maule earthquake in Chile. The cumulative afterslip distribution is obtained from the inversion of the residual surface displacements between the observed displacements from the continuous GPS (cGPS) and the ones from the forward modelling. We investigate five simulations, four with different dislocation creep parameters for the crust, slab, and upper mantle and one with elastic properties for the crust and slab, and a linear viscoelastic upper mantle for comparison. Our preferred simulation considers a weak crust since it shows the best fit to the cumulative cGPS postseismic displacements, a good fit to the time-series, and, in particular, a good spatial correlation between afterslip and aftershock activity. In this simulation, most of the viscoelastic relaxation occurs in the continental lower crust beneath the volcanic arc due to dislocation creep processes. The resulting afterslip pattern from the inversion is reduced at depths >60 km, which correlates to the low cumulative seismic moment that is released from aftershocks at these depths. Furthermore, the cumulative afterslip moment release from this simulation corresponds to 10% of the main shock in six years, which is approximately half of the moment release that results from models with an elastic crust and linear viscosity in the upper mantle. We conclude that an integrated analysis by considering power-rheology with dislocation creep processes in the continental crust and upper mantle along with aftershock activity may be used to constrain location and magnitude postseismic relaxation processes better.

Tidal Deformations of Hybrid Stars with Sharp Phase Transitions and Elastic Crusts

Abstract Gravitational wave astronomy is expected to provide independent constraints on neutron-star properties, such as their equation of state. This is possible with the measurements of binary components’ tidal deformability, which alter the point-particle gravitational waveforms of neutron-star binaries. Here, we provide a first study of the tidal deformability effects due to the elasticity/solidity of the crust (hadronic phase) in a hybrid neutron star, as well as the influence of a quark-hadronic phase density jump on tidal deformations. We employ the framework of non-radial perturbations with zero frequency and study hadronic phases presenting elastic aspects when perturbed (with the shear modulus approximately 1% of the pressure). We find that the relative tidal deformation change in a hybrid star with a perfect-fluid quark phase and a hadronic phase presenting an elastic part is never larger than about 2%–4% (with respect to a perfect-fluid counterpart). These maximum changes occur when the elastic region of a hybrid star is larger than approximately 60% of the star’s radius, which may happen when its quark phase is small and the density jump is large enough, or even when a hybrid star has an elastic mixed phase. For other cases, tidal deformation changes due to an elastic crust are negligible ( ) and, therefore, unlikely to be measured even with third generation detectors. Thus, only when the size of the elastic hadronic region of a hybrid star is over half of its radius, could the effects of elasticity have a noticeable impact on tidal deformations.

Interseismic Loading of Subduction Megathrust Drives Long‐Term Uplift in Northern Chile

AbstractLarge earthquakes are the product of elastic stress that has accumulated over decades to centuries along segments of active faults. Assuming an elastic crust, one can roughly estimate the location and rate of accumulation of elastic stress. However, this general framework does not account for inelastic, irrecoverable deformation, which results in large‐scale topography. We do not know over which part of the earthquake cycle such deformation occurs. Using InSAR and GNSS measurements, we report on a potential correlation between long‐term, inelastic vertical rate and short‐term, interseismic vertical rate in northern Chile. Approximately 4% to 8% of the geodetically derived interseismic vertical rates translate into permanent deformation, suggesting that topography of the forearc builds up during the interseismic period. This observation provides a quantitative basis for an improved understanding of the interplay between short‐term and long‐term dynamics along convergent plate boundaries.

Sensitivities of Geodetic Source Analyses to Elastic Crust Heterogeneity Constrained by Seismic Tomography for the 2017 Mw 6.5 Jiuzhaigou, China, Earthquake

ABSTRACTThe 2017 Mw 6.5 Jiuzhaigou earthquake (JE) struck a rugged area of the Jiuzhaigou Valley in eastern Tibet that has experienced frequent seismic activity over the last few decades. We use finite‐element models (FEMs) and Sentinel‐1 Interferometric Synthetic Aperture Radar observations to characterize the earthquake source. The FEM domain accommodates a heterogeneous (HET) distribution of realistic crustal materials inferred by regional seismic tomography data. The HET‐derived source configurations yield a significantly smaller misfit, at the 95% confidence level, than that estimated for a homogeneous (HOM) half‐space. The former generally requires a lower degree of smoothing constraint, highlighting that the HET solutions are systematically more compatible with the surface observations than the HOM solutions. The magnitudes of induced Coulomb failure stress change (ΔCFS) estimated by the HET solution drastically differ (by >0.1 MPa) from those calculated by the HOM solution. The postearthquake stability of near‐field faults is generally overestimated by the HOM estimations, whereas some localities of negative ΔCFSHOM are predicted with positive ΔCFSHET. These results highlight the sensitivities of both slip and stress estimations to the complexity of the adopted elastic modeling domain, leading to more accurate aftershock hazard assessments. The HET‐resolved seismic rupture reveals two major slip asperities of magnitude up to 0.83 m distributed along the fault strike, which is coherent with the aftershock distribution. Two aftershock clusters are consistently found near or below these two peak‐slip zones, which are imaged by the HET model but absent in the HOM solution. The JE hypocenter and aftershocks are bounded below by a negative velocity anomaly (ΔVP, ΔVS down to −4%) at ∼18 km depth. Such low‐velocity layers of reduced strength may be relevant to the vertical distribution of seismicity and earthquake slip, which provide insights into assessing the seismic hazards and aftershock‐prone areas of the eastern Tibetan margin.

MODERN MOVEMENTS OF THE CRUST SURFACE IN GORNY ALTAI FROM GPS DATA

In 2000–2017, the GPS technology was first applied to study inter-seismic, co-seismic and post-seismic processes in the crust of the Altai Mountains (Gorny Altai). Our study aims at investigating the fields of displacement and deformation in the Gorny Altai region as a part of Asia.The 3D displacement fields are reconstructed for the period before the M 7.3 Chuya earthquake that occurred in the southern sector of the Altai GPS network (49° to 55°N, and 81° to 89°E)on 27 September 2003.Anomalous behavior features are discovered in the displacement orientations, as well as in the distribution of velocities and deformation in the zone of the future earthquake.The spatial displacement pattern defined for the period of co-seismic displacements corresponds to the right-lateral strike-slip along the vertical fault. The fault depth is estimated using the elastic model and the experimental data (change in displacement from 0.30 m to 0.02 m at the distances of 14 km and 84 km from the fault, respectively); it amounts to 8–10 km.The co-seismic deformation field is investigated.In the post-seismic stage (2004–2017), displacements revealedin the epicentral zone show the right-lateral strike-slip along the fault at the rate of 2 mm/yr. Therefore, two-layer viscoelastic models can be considered. The estimated viscosity of the lower crust ranges from 6×1019to 3×1020Pa×s, and the elastic upper crust thickness is 25 km. Analyzed are modern movements in the Gorny Altai region outside the Chuya earthquake area.The results of our study show that modern horizontal displacements occur in the NNW direction at the rate of 1.1 mm/yr, which is twice lower than the displacement rate before the earthquake.

Lithospheric flexure and rheology determined by climate cycle markers in the Corinth Rift

Geomorphic strain markers accumulating the effects of many earthquake cycles help to constrain the mechanical behaviour of continental rift systems as well as the related seismic hazards. In the Corinth Rift (Greece), the unique record of onshore and offshore markers of Pleistocene ~100-ka climate cycles provides an outstanding possibility to constrain rift mechanics over a range of timescales. Here we use high-resolution topography to analyse the 3D geometry of a sequence of Pleistocene emerged marine terraces associated with flexural rift-flank uplift. We integrate this onshore dataset with offshore seismic data to provide a synoptic view of the flexural deformation across the rift. This allows us to derive an average slip rate of 4.5–9.0 mm·yr−1 on the master fault over the past ~610 ka and an uplift/subsidence ratio of 1:1.1–2.4. We reproduce the observed flexure patterns, using 3 and 5-layered lithospheric scale finite element models. Modelling results imply that the observed elastic flexure is produced by coseismic slip along 40–60° planar normal faults in the elastic upper crust, followed by postseismic viscous relaxation occurring within the basal lower crust or upper mantle. We suggest that such a mechanism may typify rapid localised extension of continental lithosphere.

Magma chambers: what we can, and cannot, learn from volcano geodesy.

Geodetic observations on volcanoes can reveal important aspects of crustal magma chambers. The rate of decay of deformation with distance reflects the centroid depth of the chamber. The amplitude of the deformation is proportional to the product of the pressure change and volume of the reservoir. The ratio of horizontal to vertical displacement is sensitive to chamber shape: sills are efficient at generating vertical displacement, while stocks produce more horizontal deformation. Geodesy alone cannot constrain important parameters such as chamber volume or pressure; furthermore, kinematic models have no predictive power. Elastic response combined with influx proportional to pressure gradient predicts an exponentially decaying flux, leading to saw-tooth inflation cycles observed at some volcanoes. Yet many magmatic systems exhibit more complex temporal behaviour. Wall rock adjacent to magma reservoirs cannot behave fully elastically. Modern conceptual models of magma chambers also include cumulate and/or mush zones, with potentially multi-level melt lenses. A viscoelastic shell surrounding a spherical magma chamber significantly modifies the predicted time-dependent response; post-eruptive inflation can occur without recharge if the magma is sufficiently incompressible relative to the surrounding crust (Segall P. 2016 J. Geophys. Res. Solid Earth, 121, 8501-8522). Numerical calculations confirm this behaviour for both oblate and prolate ellipsoidal chambers surrounded by viscoelastic aureoles. Interestingly, the response to a nearly instantaneous pressure drop during an explosive eruption can be non-monotonic as the rock around the chamber relaxes at different rates. Pressure-dependent recharge of a non-Newtonian magma in an elastic crust leads to an initially high rate of inflation which slows over time; behaviour that has been observed in some magmatic systems. I close by discussing future challenges in volcano geodesy. This article is part of the Theo Murphy meeting issue 'Magma reservoir architecture and dynamics'.

Incompressible analytical models for spinning-down pulsars

Abstract We study a class of Newtonian models for the deformations of non-magnetised neutron stars during their spin-down. All the models have an analytical solution which allows to easily grasp the dependence of the strain on the star’s main physical quantities, such as radius, mass, and crust thickness. We first use the model proposed by Franco, Link, and Epstein that depicts the star as made of a fluid core and an elastic crust with the same density, to compare the response to a decreasing centrifugal force on stars having different masses and equations of state. We find that the strain angle is peaked at the equator and its maximum value decreases as a function of the mass. Afterwards, we introduce a second, more refined, model in which the core and the crust have different densities, and the gravitational potential of the deformed body is self-consistently accounted for. The strain angle is still a decreasing function of the stellar mass, but now its maximum value is typically peaked at the poles and is larger (by a factor of four) than the corresponding value in the one-density model. Finally, within the present analytic approach, we evaluate the impact of the Cowling approximation: when the perturbations of the gravitational potential are neglected, we find an underestimation of the centrifugal effect on the star, since the strain angle is about 40% of the one obtained with the complete model.

On the Mechanical Effects of Poroelastic Crystal Mush in Classical Magma Chamber Models

AbstractImproved constraints on the mechanical behavior of magma chambers is essential for understanding volcanic processes; however, the role of crystal mush on the mechanical evolution of magma chambers has not yet been systematically studied. Existing magma chamber models typically consider magma chambers to be isolated melt bodies surrounded by elastic crust. In this study, we develop a physical model to account for the presence and properties of crystal mush in magma chambers and investigate its impact on the mechanical processes during and after injection of new magma. Our model assumes the magma chamber to be a spherical body consisting of a liquid core of fluid magma within a shell of crystal mush that behaves primarily as a poroelastic material. We investigate the characteristics of time‐dependent evolution in the magma chamber, both during and after the injection, and find that quantities such as overpressure and tensile stress continue to evolve after the injection has stopped, a feature that is absent in elastic (mushless) models. The time scales relevant to the postinjection evolution vary from hours to thousands of years, depending on the micromechanical properties of the mush, the viscosity of magma, and chamber size. We compare our poroelastic results to the behavior of a magma chamber with an effectively viscoelastic shell and find that only the poroelastic model displays a time scale dependence on the size of the chamber for any fixed mush volume fraction. This study demonstrates that crystal mush can significantly influence the mechanical behaviors of crustal magmatic reservoirs.

Episodic kinematics in continental rifts modulated by changes in mantle melt fraction.

Oceanic crust is created by the extraction of molten rock from underlying mantle at the seafloor 'spreading centres' found between diverging tectonic plates. Modelling studies have suggested that mantle melting can occur through decompression as the mantle flows upwards beneath spreading centres, but direct observation of this process is difficult beneath the oceans. Continental rifts, however-which are also associated with mantle melt production-are amenable to detailed measurements of their short-term kinematics using geodetic techniques. Here we show that such data can provide evidence for an upwelling mantle flow, as well as information on the dimensions and timescale of mantle melting. For North Island, New Zealand, around ten years of campaign and continuous GPS measurements in the continental rift system known as the Taupo volcanic zone reveal that it is extending at a rate of 6-15 millimetres per year. However, a roughly 70-kilometre-long segment of the rift axis is associated with strong horizontal contraction and rapid subsidence, and is flanked by regions of extension and uplift. These features fit a simple model that involves flexure of an elastic upper crust, which is pulled downwards or pushed upwards along the rift axis by a driving force located at a depth greater than 15 kilometres. We propose that flexure is caused by melt-induced episodic changes in the vertical flow forces that are generated by upwelling mantle beneath the rift axis, triggering a transient lower-crustal flow. A drop in the melt fraction owing to melt extraction raises the mantle flow viscosity and drives subsidence, whereas melt accumulation reduces viscosity and allows uplift-processes that are also likely to occur in oceanic spreading centres.

A possible origin of the Galactic Center magnetar SGR 1745–2900

Since there is a large population of massive O/B stars and putative neutron stars (NSs) located in the vicinity of the Galactic Center (GC), intermediate-mass X-ray binaries (IMXBs) constituted by an NS and a B-type star probably exist there. We investigate the evolutions of accreting NSs in IMXBs (similar to M82 X-2) with a companion and orbital period d. By adopting a mildly super-Eddington rate for the early Case B Roche-lobe overflow (RLOF) accretion, we find that only in accreting NSs with quite elastic crusts (slippage factor s = 0.05) can the toroidal magnetic fields be amplified within 1 Myr, which is assumed to be the longest duration of the RLOF. These IMXBs will evolve into NS+white dwarf (WD) binaries if they are dynamically stable. However, before the formation of NS+WD binaries, the high stellar density in the GC will probably lead to frequent encounters between the NS+evolved star binaries (in post-early Case B mass transfer phase) and NSs or exchange encounters with other stars, which may produce single NSs. These NSs will evolve into magnetars when the amplified poloidal magnetic fields diffuse out to the NS surfaces. Consequently, our results provide a possible explanation for the origin of the GC magnetar SGR 1745–2900. Moreover, the accreting NSs with will evolve into millisecond pulsars (MSPs). Therefore, our model reveals that the GC magnetars and MSPs could both originate from a special kind of IMXB.

Accumulation of elastic strains in the upper crust on the locked transform faults and the tectonomagnetic effect

Plate tectonics only allows small deformations in the lithospheric plates. The laboratory experiments with the rock specimens show that the creep is transient when the creep strain is at most 1%. Hence, if we assume that the creep strain in the lithospheric plates is below this threshold, the creep is transient. The present paper addresses the role of the elastic, brittle (pseudo-plastic), and creep rheology of the lithosphere during the accumulation of elastic shear strains on the locked faults in the Earth’s crust, i.e., during the process of preparation of the earthquakes. The effective viscosity characterizing the transient creep is lower than that under the steady-state creep and it depends on the characteristic time of a given process. The characteristic duration of the stress and strain accumulation process in the vicinity of the locked faults is a few dozen years. On these time intervals, the thin upper crustal layer behaves as brittle; the underlying layer behaves as elastic (it is just this layer which accommodates stress accumulation leading to the earthquake), whereas the transient creep is predominant in the lower crust and mantle lithosphere. Transient creep entails nonlinear time dependence of the strains arising in the vicinity of the locked fault in the elastic crust. The perturbations in the magnetic field induced by these strains can be treated as the magnetic precursor of the earthquake.

Geodetically constrained models of viscoelastic stress transfer and earthquake triggering along the North Anatolian fault

AbstractOver the past 80 years, 8 MW > 6.7 strike‐slip earthquakes west of 40° longitude have ruptured the North Anatolian fault (NAF) from east to west. The series began with the 1939 Erzincan earthquake in eastern Turkey, and the most recent 1999 MW = 7.4 Izmit earthquake extended the pattern of ruptures into the Sea of Marmara in western Turkey. The mean time between seismic events in this westward progression is 8.5 ± 11 years (67% confidence interval), much greater than the timescale of seismic wave propagation (seconds to minutes). The delayed triggering of these earthquakes may be explained by the propagation of earthquake‐generated diffusive viscoelastic fronts within the upper mantle that slowly increase the Coulomb failure stress change ( CFS) at adjacent hypocenters. Here we develop three‐dimensional stress transfer models with an elastic upper crust coupled to a viscoelastic Burgers rheology mantle. Both the Maxwell (ηM = 4 × 1018−1 × 1019 Pa s) and Kelvin (ηK = 1 × 1018−1 × 1019 Pa s) viscosities are constrained by studies of geodetic observations before and after the 1999 Izmit earthquake. We combine this geodetically constrained rheological model with the observed sequence of large earthquakes since 1939 to calculate the time evolution of CFS changes along the North Anatolian fault due to viscoelastic stress transfer. Apparent threshold values of mean CFS at which the earthquakes in the eight decade sequence occur are between ∼0.02 to ∼3.15 MPa and may exceed the magnitude of static CFS values by as much as 177%. By 2023, we infer that the mean time‐dependent stress change along the northern NAF strand in the Marmara Sea near Istanbul, which may have previously ruptured in 1766, may reach the mean apparent time‐dependent stress thresholds of the previous NAF earthquakes.