# Time-dependent geoid anomalies at subduction zones due to the seismic cycle

Time-dependent geoid anomalies at subduction zones due to the seismic cycle Abstract We model the geoid anomalies excited during a megathrust earthquake cycle at subduction zones, including the interseismic phase and the contribution from the infinite series of previous earthquakes, within the frame of self-gravitating, spherically symmetric, compressible, viscoelastic Earth models. The fault cuts the whole 50 km lithosphere, dips 20°, and the slip amplitude, together with the length of the fault, are chosen in order to simulate an Mw = 9.0 earthquake, while the viscosity of the 170 km thick asthenosphere ranges from 1017 to 1020 Pa s. On the basis of a new analysis from the Correspondence Principle, we show that the geoid anomaly is characterized by a periodic anomaly due to the elastic and viscous contribution from past earthquakes and to the back-slip of the interseismic phase, and by a smaller static contribution from the steady-state response to the previous infinite earthquake cycles. For asthenospheric viscosities from 1017–1018 to 1019–1020 Pa s, the characteristic relaxation times of the Earth model change from shorter to longer timescales compared to the 400 yr earthquake recurrence time, which dampen the geoid anomaly for the higher asthenospheric viscosities, since the slower relaxation cannot contribute its whole strength within the interseismic cycle. The geoid anomaly pattern is characterized by a global, time-dependent positive upwarping of the geoid topography, involving the whole hanging wall and partially the footwall compared to the sharper elastic contribution, attaining, for a moment magnitude Mw = 9.0, amplitudes as high as 6.6 cm for the lowermost asthenospheric viscosities during the viscoelastic response compared to the elastic maximum of 3.8 cm. The geoid anomaly vanishes due to the back-slip of the interseismic phase, leading to its disappearance at the end of the cycle before the next earthquake. Our results are of importance for understanding the post-seismic and interseismic geoid patterns at subduction zones. Seismic cycle, Subduction zone processes, Dynamics: gravity and tectonics, Dynamics: seismoctectonics 1 INTRODUCTION The 2004 Sumatra, 2010 Maule and 2011 Tohoku megathrust earthquakes have caused the flourishing of studies on the gravity signatures caused by the earthquake-induced mass redistribution, nowadays detectable by GRACE and GOCE satellite missions. The co-seismic gravity anomalies caused by earthquakes with moment magnitude greater than 8.5 have been examined by means of long-wavelength gravity data from GRACE (Gross & Chao 2001; Sun & Okubo 2004; Matsuo & Heki 2011) and explained in terms of mass redistribution, volume variations and topography perturbations by De Linage et al. (2009), Cambiotti et al. (2011a) and Broerse et al. (2011). For the case of the 2011 Tohoku earthquake, GRACE data have been used for estimating its focal mechanism and hypocentre (Wang et al. 2012; Cambiotti & Sabadini 2013), and, jointly with GPS and GOCE data, to improve the constraint of its slip distribution over fault surface (Fuchs et al. 2016). In addition to the co-seismic signature, the post-seismic response associated to stress relaxation by viscous flow of the 2004 Sumatra earthquake has been investigated in a series of work (Pollitz et al. 2006; Chen et al. 2007; Han et al. 2008; Broerse et al. 2015; Tanaka et al. 2015) with the aim of constraining the viscosity stratification of the lithosphere and of the asthenosphere. Due to the monthly resolution of the GRACE satellite mission, estimates of the co-seismic jump at the earthquake time in the data time series are usually performed jointly with the estimate of the post-seismic response after the earthquake. Nevertheless, the interpretation of the latter requires a full understanding of all the contributions from the seismic cycle. Indeed, the time-dependent signature after the earthquake, in addition to the post-seismic response to the earthquake, is affected also by stress relaxation by viscous flow from past seismicity, as well as by strain accumulation due to the locking of the fault during the interseismic phase. As it concerns the topography changes at subduction zones, this issue has been first addressed by Savage (1983), accounting for both strain accumulation and release within the frame of viscoelastic, flat Earth models (Thatcher & Rundle 1979). According to Savage & Burford (1973) and Savage (1983) and as shown by Matsu’ura & Sato (1989) and Sato & Matsu’ura (1993) for old subduction zones, the long term average motion between the overriding and subducting plates is responsible for static vertical displacements, as for the peculiar bipolar gravity, geoid and topography patterns at subduction zones (Zhong & Gurnis 1992; Song & Simons 2003; Marotta et al. 2006), but it does not contribute to observable changes through time. In addition to the instantaneous elastic response at the earthquake time, changes in time during the interseismic phase are the results of the post-seismic response to the last earthquake, as well as to all the past seismicity and to the viscoelastic response to strain accumulation. Only for young subduction zones (i.e. for timescales smaller than the characteristic relaxation time of the high viscosity lithosphere) this relative motion yields a growth rate of the topography of island arc-trench systems. In this work, we will consider geoid anomalies at subduction zones due to both seismic and interseismic processes. Particularly, in order to simulate the geoid anomalies due to the locking of the fault during the interseismic phase, we follow the approach first proposed by Savage & Burford (1973) for strike-slip faults and, then, generalized to dip-slip faults in Savage (1983) (see also Segall 2010). According to this approach, the strain accumulation due to the locking of the fault is taken into account by imposing a constant in time slip rate on the fault system associated to the plate boundary that is opposite to the average slip rate necessary to accommodate the relative motion between plates. Compared to Savage (1983), we implement a compressible viscoelastic Earth model which is spherically symmetric and self-gravitating, thus allowing for a self-consistent treatment of the geoid anomalies, and builds on the elastic and density stratifications of the preliminary reference Earth model (PREM, Dziewonski & Anderson 1981). Furthermore, rather than considering separately the response to each earthquake and to the strain accumulation, we take advantage of the Correspondence Principle between the viscoelastic problem in the time domain and the equivalent elastic problem in the Laplace domain (Tanaka et al. 2006; Cambiotti et al. 2011c) in order to evaluate analytically the static geoid anomaly over which the periodic changes during the interseismic phase are superimposed. 2 SLIP TIME HISTORY As shown in Fig. 1, we approximate the fault surface at the subduction zone with a planar fault dipping with an angle α = 20° and we assume that the locked portion during the interseismic phase extends through the whole 50 km thick outermost elastic layer of the Earth model. Figure 1. View largeDownload slide Schematicdepicting the 50 km thick outermost elastic layer, the 170 km thick asthenosphere and upper and lower mantle below. The locked portion of the fault, indicated by the dashed line, dips with angle α = 20° and extends from the Earth surface down to a depth Z = 50 km. Its width is W = Z/sin α = 146 km. The coordinate x perpendicular to the line of strike, that is the distance from the trench over the Earth surface, is also shown. Positive and negative values identify points on the sides of the subducting (foot wall) and (hanging wall) plates, as indicated by the dash-dotted arrow. Figure 1. View largeDownload slide Schematicdepicting the 50 km thick outermost elastic layer, the 170 km thick asthenosphere and upper and lower mantle below. The locked portion of the fault, indicated by the dashed line, dips with angle α = 20° and extends from the Earth surface down to a depth Z = 50 km. Its width is W = Z/sin α = 146 km. The coordinate x perpendicular to the line of strike, that is the distance from the trench over the Earth surface, is also shown. Positive and negative values identify points on the sides of the subducting (foot wall) and (hanging wall) plates, as indicated by the dash-dotted arrow. For the sake of simplicity and with the aim of providing a first estimate of the geoid anomaly due to seismic and interseismic processes at subduction zones, we assume that thrust earthquakes occur at a constant recurrence period, T, and break uniformly the locked portion of the fault, with the same amount U of dip-slip at each earthquake. We thus write the stick-slip time history of the earthquake sequence as follows:   \begin{eqnarray} \delta u_e(t) = U\, \sum _{n=1}^{\infty } H(t-n\, T), \end{eqnarray} (1)where t is the time and H is the Heaviside function. Herein we assume that the first earthquake occurs at time t = T and, by convention, we start to count the seismic cycles from this time, This means that the nth seismic cycle is from n T to (n + 1) T. The initial period, from 0 to T, must be regarded as an arbitrary loading phase before the occurrence of the first earthquake. In addition to the stick-slip from each earthquake and according to Savage (1983), we consider also a back-slip distributed uniformly over the locked portion of the fault which grows linearly in time, starting from t = 0,   \begin{eqnarray} \delta u_b(t) = - \chi \, v\, t\, H(t) \end{eqnarray} (2)where χ is the seismic coupling coefficient and v is the convergence velocity between the subducting and overriding plates. This back-slip accounts for the strain accumulation process due to plate tectonics and constitutes the supplemental forcing that, once superimposed to the contribution from the steady state subduction (i.e. the long-term average motion that is not responsible for observable changes), describes the locking of the fault during the interseismic phase. Depending on the seismic coupling coefficient, this locking can be complete (χ = 1) or partial (χ < 1). For consistency, we equal the back-slip accumulated during one period, δub(t + T) − δub(t) = −χ v T to the stick-slip of one earthquake changed in sign, −U, and we obtain the following relation:   \begin{eqnarray} U = \chi \, v\, T. \end{eqnarray} (3) The dip-slip time history for which we will compute the geoid anomalies is the sum of the stick-slip and back-slip time histories:   \begin{eqnarray} \delta u(t) = \delta u_e(t)+ \delta u_b(t) = U\, f(t) \end{eqnarray} (4)with   \begin{eqnarray} f(t) = -\frac{t}{T}+\sum _{n=1}^{\infty } H(t-n\, T). \end{eqnarray} (5)In view of the arbitrariness of the choice of the initial time, the simulations that we show are significant only in the limit for t → ∞, after an infinite number of seismic cycles. We note that eq. (5) is equivalent to eq. (3) of Matsu’ura & Sato (1989) and, as in our case, it describes the slip history responsible for the periodic motion associated to the infinite sequence of seismic cycles within the seismogenic zone. Matsu’ura & Sato (1989) further account for the relative motion across the whole plate boundary between the subducting and overriding plates, their eq. (2). Nevertheless, as discussed by Matsu’ura & Sato (1989) and Sato & Matsu’ura (1993) for the case of old subduction zones, this relative motion is responsible only for a static contribution and, so, it does not yield any observable changes, in agreement with the assumptions underlying the back-slip approach proposed by Savage (1983). In this work, in order to focus on the fastest periodic changes within one seismic cycle, we do not consider the case of young subduction zones where this relative motion yields a growth rate of the topography of island arc-trench systems from one seismic cycle to the next one. 3 VISCOELASTIC RESPONSE From eq. (A6) and as discussed in Appendix A, within the framework of self-gravitating spherically symmetric viscoelastic Earth models with Maxwell rheology (Tanaka et al. 2006; Cambiotti et al. 2011c; Cambiotti & Sabadini 2015), the geoid anomaly at a point of the Earth surface, G, can be expressed as follows:   \begin{eqnarray} G(t) = (K\star \delta u)(t), \end{eqnarray} (6)where ⋆ stands for the time convolution and K is the viscoelastic Green function for a unitary, impulsive dip-slip over the whole locked portion of the fault. In order to compute the viscoelastic Green function, we consider a spherically symmetric Earth model based on PREM (table 3 of Dziewonski & Anderson (1981)) characterized by a 50 km thick elastic outermost layer and a 170 km thick (up to 220 km depth) asthenosphere, as portrayed in Fig. 1, with viscosity νa varied from 1017 to 1020 Pa s. The viscosities of the upper (from 220 km depth to the 670 km discontinuity) and lower (from the 670 km discontinuity to the core mantle boundary) mantle are set to 1021 and 3 × 1022 Pa s. This rheological stratification is consistent with most glacial isostatic adjustment and mantle convection studies (Cambiotti et al. 2011b; Tosi et al. 2005) and includes a low viscosity asthenosphere that mainly controls the relaxation of the earthquake-induced stress by viscous flow (Broerse et al. 2015). Furthermore, the initial density and bulk modulus are slightly modified in order to satisfy the Williamson–Adams equation (Wolf 1991; Wolf & Kaufmann 2000) and avoid Rayleigh–Taylor instabilities or compositional modes at large timescales (Plat & Jüttner 1995; Cambiotti & Sabadini 2010; Cambiotti et al. 2013). As shown by eqs (A12), (B10), (B11), (B13) and (A13a) and as discussed in Appendices A and B, eq. (6) with the dip-slip time history given by eq. (4) can be decomposed into three contributions:   \begin{eqnarray} G(t) = U\, \big (K_E\, f(t)+P(t)+C(t)\big ), \end{eqnarray} (7)where KE is the elastic Green function, P(t) is a periodic function with period T  \begin{eqnarray} P(t) = P(t+n\, T) \end{eqnarray} (8)such that P(n T) = 0, and C(t) has the following asymptotic limit   \begin{eqnarray} \lim _{t\rightarrow \infty }C(t) = C_\infty . \end{eqnarray} (9)The periodic function P(t) and this asymptotic limit can be expressed analytically, eqs (B14) and (B11), once the Laplace transform of the viscoelastic Green function has been evaluated at specific values of the Laplace variable. In view of the periodicity in time of the functions f(t) and P(t), the response of the Earth during the nth seismic cycle can be rewritten as   \begin{eqnarray} G_n(t) &= G(t+n\, T) = U \left(E(t)+P(t)+S_n(t)\right) & \end{eqnarray} (10)with t ∈ [0, T) and   \begin{eqnarray} E(t) = \left(1-t/T\right)\, K_E \end{eqnarray} (11a)  \begin{eqnarray} S_n(t) = C(t+n\, T)-K_E. & \end{eqnarray} (11b) We note that Sn(t) depends on the seismic cycle, while E(t) and P(t) do not. On the other hand, in view of eq. (9), we have   \begin{eqnarray} S_\infty = \lim _{n\rightarrow \infty } S_n(t) = C_\infty -K_E & \end{eqnarray} (12)and, then, the geoid anomaly after an infinite number of seismic cycles becomes a fully periodic function   \begin{eqnarray} G_\infty (t) = \lim _{n\rightarrow \infty }G_n(t)=U\, \big (E(t)+P(t)+S_\infty \big ). & \end{eqnarray} (13)From the latter formulation, we realize that E(t) and P(t) control the periodic changes of the geoid anomaly during one seismic cycle due to the instantaneous elastic response of the Earth and to the stress relaxation by viscous flow. Particularly, these periodic changes account for the instantaneous elastic response at the earthquake time, t = 0, and for the following changes during the interseismic phase which become zero just before the next earthquake, t = T−,   \begin{eqnarray} E(0)+P(0) = K_E \end{eqnarray} (14a)  \begin{eqnarray} E(T^-)+P(T^-) = 0. \end{eqnarray} (14b)We note also that the instantaneous elastic response, E(t), starts with the elastic response to stick-slip at the earthquake time and decreases linearly in time to zero just before the next earthquake. The constant S∞, instead, controls the static geoid anomaly over which these periodic changes are superimposed. This contribution results from stress relaxation by viscous flow due to the previous (infinite) seismic cycles and corresponds to the geoid anomaly evaluated just before the next earthquake that, as implied by eq. (12), reaches an asymptotic value after an infinite number of seismic cycles   \begin{eqnarray} G_\infty (T^{-}) = U\, S_\infty . \end{eqnarray} (15)This static geoid anomaly must not be confused with that caused by the steady state subduction. Our model, in fact, only focus on the strain accumulation and release process and does not account for the long-term relative motion between the subducting and overriding plates as in Matsu’ura & Sato (1989) and Sato & Matsu’ura (1993) and for the density heterogeneities of the subducting slab into the mantle and the induced mantle flow as in Ricard et al. (1993). Rather, as we are going to show in Section 5, it is a much smaller contribution resulting from the specific setting of the dip-slip time history describing the strain accumulation and release and that, although small and negligible, should be added to the static geoid anomaly caused by the steady state subduction. This analysis of the geoid anomaly after an infinite number of seismic cycles has allowed us to understand the physical meaning of E(t), P(t) and S∞. This understanding can now be used to gain insight into the response of the Earth after a finite number of seismic cycles. In this case, there is a net gain of the geoid anomaly accumulated during one interseismic period   \begin{eqnarray} G_{n+1}(t)-G_n(t) = S_{n+1}(t)-S_n(t) \end{eqnarray} (16)which involves only the function Sn(t). This net gain results from the fact that the stress relaxation by viscous flow due to all previous (finite) earthquakes and to the back-slip simulating the locking of the fault during the interseismic phase has not reached the steady state yet. Increasing the number of seismic cycles, in the limit for n → ∞, this net gain tends to zero and, then, there are no longer differences between a seismic cycle and the next one. 4 THE CASE OF A SIMPLE RELAXATION SPECTRUM In order to gain insight into the interaction between the stress relaxation and the dip-slip time history of the sequence of seismic cycles, it is useful to consider a heuristic Maxwell Earth model with a discrete relaxation spectrum. Within this assumption, the viscoelastic Green function can be expressed as follows:   \begin{eqnarray} K(t) = K_E\, \delta (t)+H(t)\, \sum _j k_j\, e^{-t/T_j} \end{eqnarray} (17)where δ(t) is the Dirac delta function, and kj and Tj are the residue and characteristic relaxation time of the jth relaxation mode. This specific form allows us to obtain the analytical expression for the time-dependent geoid anomaly. Particularly, as discussed in Appendix C, eq. (C5a), the periodic change due only to stress relaxation by viscous flow reads   \begin{eqnarray} P(t) = \sum _j r_j\, p_j(t) \end{eqnarray} (18)where rj = Tj Kj is the response to a stick-slip earthquake in the fluid limit, for t → ∞, associated to the jth relaxation mode (the so called strength) and pj(t) is a non-dimensional function describing the time evolution of this contribution during the seismic cycle   \begin{eqnarray} p_j(t) = \frac{1-e^{-t/T_j}}{1-e^{-T/T_j}}-\frac{t}{T}. \end{eqnarray} (19)This function is zero at the earthquake time, t = 0, increases to a maximum value $$p_j^*$$ at $$t=t^*_j$$, and then decreases to zero at the end of the seismic cycle, t = T. It is straightforward to obtain   \begin{eqnarray} p_j^* = \frac{1}{1-e^{-T/T_j}} +\frac{t_j^*-T_j}{T} \end{eqnarray} (20a)  \begin{eqnarray} t_j^* = T_j\, \left[\log \left(\frac{T}{T_j}\right)-\log \left(1-e^{-T/T_j}\right)\right]. \end{eqnarray} (20b) Fig. 2(a) shows the evolution in time of the function pj(t) for different relaxation times Tj = 0.01, 0.03, 0.1, 0.3, 1 and 3 T. We note that, increasing the relaxation times, the maximum amplitude $$p^*_j$$ decreases, while the time $$t_j^*$$ at which this maximum amplitude is attained increases. Fig. 2(b) shows also the maximum amplitude $$p_j^*$$ and the time $$t_j^*$$ at which this value is attained as function of the relaxation time Tj. We note that   \begin{eqnarray} \lim _{T_j/T\rightarrow 0} p_j^* = 1 \end{eqnarray} (21a)  \begin{eqnarray} \lim _{T_j/T\rightarrow \infty } p_j^* = 0. \end{eqnarray} (21b) This means that, for relaxation times smaller than the interseismic period, Tj ≪ T, the contribution of the relaxation mode corresponds almost to its strength, while for greater relaxation times, Tj ≫ T, its contribution is almost filtered out. This damping is already effective for relaxation times comparable with the interseismic period and, particularly, $$p_j^*\approx 1/2$$ for Tj ≈ T/5. As it concerns the time $$t_j^*$$ at which the contribution attains its maximum, it is always smaller than T/2 and this upper bound is reached only for relaxation times greater than the interseismic period, in the limit for Tj/T → ∞. In contrast, for relaxation time smaller than the interseismic period, Tj ≪ T, the time $$t_j^*$$ is about the relaxation time Tj scaled by log (T/Tj). Figure 2. View largeDownload slide (a) The contribution of the jth relaxation mode per unit strength during one seismic cycles of period, pj(t), for different relaxation times Tj = 0.01, 0.1, 0.3, 1 and 3 T. (b) The maximum amplitude of this contribution (left vertical axis, solid line) and the time at which this maximum is reached (right vertical axis, dashed line) as function of the relaxation time. Figure 2. View largeDownload slide (a) The contribution of the jth relaxation mode per unit strength during one seismic cycles of period, pj(t), for different relaxation times Tj = 0.01, 0.1, 0.3, 1 and 3 T. (b) The maximum amplitude of this contribution (left vertical axis, solid line) and the time at which this maximum is reached (right vertical axis, dashed line) as function of the relaxation time. 5 RESULTS We consider an interseismic period of T = 400 yr and, from eq. (3), assuming a partial seismic coupling of χ = 0.5 and a convergence velocity of v = 10 cm/yr−1, we set the stick-slip of each earthquake at U = 20 m. Considering a fault segment of (along strike) length L = 260 km and the shear modulus profile of PREM, each earthquake has moment magnitude Mw = 9.0. In the following we will focus first on the periodic geoid anomaly during one seismic cycle due to stress relaxation by viscous flow, U P(t), and we will consider the results for different asthenospheric viscosities in order to quantify the impact of the latter on the time-dependent pattern of this contribution. Then, once we will have clarified this issue, we will consider the whole periodic geoid anomaly, U (E(t) + P(t)), including the instantaneous elastic response, U E(t). Fig. 3 shows the periodic geoid anomaly during one seismic cycle due to stress relaxation by viscous flow, U P(t), along a transect over the Earth surface centred with the fault and perpendicular to the line of strike, according to Fig 1. This simulation is run for four asthenospheric viscosities (νa = 1017, 1018, 1019 and 1020 Pa s) and the geoid anomaly is shown as function of the distance from the trench, x, as depicted in Fig. 1. This periodic contribution is zero at the earthquake time, t = 0, and manifests itself only at later times according to the characteristic relaxation times of the Maxwell Earth model. For νa = 1017 Pa s, in the early phase of the seismic cycle, we have a positive geoid anomaly centred at around x = −150 km on the side of the overriding plate. This positive anomaly extends for about 600 ∼ 1000 km perpendicularly to the line of strike, mainly over the overriding plate. Elsewhere, instead, the geoid anomaly is negative, although small, less than a few millimetres at all times. The maximum geoid anomaly is about 1.7 cm already at t = 0.1 yr (about one month after the earthquake) and grows in time until t = 1 yr, when there is a maximum geoid anomaly of about 5.2 cm. At later times, the geoid anomaly begins to decrease and, at the same time, a local (positive) minimum close to the trench, at around x = −50 km on the side of the overriding plate, appears between two (positive) maxima, one on the side of the overriding plate, at around x = −200 km, and the other on the side of the subducting plate, at around x = 150 km. The main geoid anomaly remains positive in proximity of the trench until t = 10 yr and, successively, becomes negative in correspondence of the local minimum, reaching a minimum value of about −1.8 mm at y = 70 yr. Then, the amplitude of the geoid anomaly continues to decrease for all the distances from the trench until it is zero everywhere at the end of the seismic cycle. Figure 3. View largeDownload slide Periodic geoid anomaly during one seismic cycle due to stress relaxation by viscous flow, U P(t), along a transect over the Earth surface centred with the fault and perpendicular to line of strike. The interseismic period and the stick-slip of each earthquake are set to T = 400 yr and U = 20 m. The locked fault crosses the whole elastic layer, down to a depth Z = 50 km. Each panel shows the geoid anomaly as function of the distance from the trench, x, at increasing time (from top to bottom and from left to right; see the label in each panel) for asthenosphere viscosities of 1017, 1018, 1019 and 1020 Pa s (solid, dashed, dashed-dotted and dotted lines, respectively). The vertical grey line indicates the trench, at x = 0. Figure 3. View largeDownload slide Periodic geoid anomaly during one seismic cycle due to stress relaxation by viscous flow, U P(t), along a transect over the Earth surface centred with the fault and perpendicular to line of strike. The interseismic period and the stick-slip of each earthquake are set to T = 400 yr and U = 20 m. The locked fault crosses the whole elastic layer, down to a depth Z = 50 km. Each panel shows the geoid anomaly as function of the distance from the trench, x, at increasing time (from top to bottom and from left to right; see the label in each panel) for asthenosphere viscosities of 1017, 1018, 1019 and 1020 Pa s (solid, dashed, dashed-dotted and dotted lines, respectively). The vertical grey line indicates the trench, at x = 0. The pattern of the geoid anomaly for νa = 1018 Pa s is similar, although delayed in time by about a factor of 10 with respect to the case of νa = 1017 Pa s. This means that the relaxation times of the Maxwell Earth model are mainly controlled by the asthenospheric viscosity. The maximum geoid anomaly is about 1.7 cm after one year, t = 1 yr, and grows up to 5.2 cm at 10 yr. Afterwards, the local positive minimum close to the trench and between the two local maxima begins to appear. Then, it becomes negative at 150 yr and reaches its minimum value of −0.3 cm at 250 yr, when the amplitude of the geoid anomaly is already decreasing for all distance from the trench until it is zero everywhere at the end of the seismic cycle. The cases of νa = 1019 and 1020 Pa s are simpler because the geoid anomalies consist mainly of positive geoid anomalies in proximity of the trench, the amplitudes of which grow in time until t = 70 and 150 yr, when they reach the maximum values of 4.2 and 0.9 cm, respectively. At later times, the amplitudes of these positive geoid anomalies begin to decrease for all distances from the trench and are zero at the end of the seismic cycle. It is noteworthy that the maximum amplitudes of the geoid anomaly over the whole seismic cycle for both cases of νa = 1017 and 1018 Pa s are similar (about 5.2 cm), although they are attained at different times (t = 1 and 10 yr, respectively). In contrast, the maximum amplitude of the geoid anomaly obtained for νa = 1019 Pa s is less by about the 20 and 80 per cent than this value, only 4.2 and 0.9 cm at t = 70 and 150 yr, respectively. This reduction of the amplitude of the geoid anomaly by increasing the asthenospheric viscosity is consistent with the discussion of Section 4 and is explained considering that the main relaxation modes of the Maxwell Earth models with νa = 1017 and 1018 Pa s have relaxation modes with characteristic relaxation times smaller than the interseismic period T = 400 yr and, then, they contribute with almost their complete strengths during the interseismic phase. In contrast, the main relaxation modes of the Maxwell Earth model with νa = 1019 have larger characteristic relaxation times and their contribution is partially filtered out. This is further confirmed by the simulation obtained for asthenospheric viscosity of 1020 Pa s. In this case, the temporal evolution is similar to that obtained for νa = 1019 Pa s, but the maximum amplitude is even smaller. Having discussed the periodic geoid anomaly due to stress relaxation by viscous flow, we now consider the whole periodic geoid anomaly, U (E(t) + P(t)), including the instantaneous elastic response as shown in Fig. 4. The geoid anomaly at the earthquake time, t = 0, is the same for all asthenospheric viscosities and corresponds to the elastic response to the stick-slip earthquake, U KE. It is characterized by a positive geoid anomaly in proximity of the trench, with a maximum of about 3.8 cm on the side of the overriding plate, at around x = −40 km, and a negative geoid anomaly on the side of the overriding plate, with a minimum of about −1.3 cm at around x = −180 km. According to eq. (11a), this elastic response decreases linearly in time during the interseismic phase and yields zero at the end of the seismic cycle. Figure 4. View largeDownload slide As Fig. 3, but including the instantaneous elastic response, U (E(t) + P(t)). Figure 4. View largeDownload slide As Fig. 3, but including the instantaneous elastic response, U (E(t) + P(t)). At the early stage, the stress relaxation increases the positive geoid anomaly caused by the elastic response in proximity of the trench, and, being wider in space, compensates the negative geoid anomaly at around x = −150 km, on the side of the overriding plate. In particular, this compensation of the negative geoid anomaly due to the instantaneous elastic response occurs at t = 0.1, 1, 10 and 150 yr for νa = 1017, 1018, 1019 and 1020 Pa s. The maximum geoid anomaly attained during the seismic cycle is 6.6 cm at t = 1 yr, 6.5 cm at 10 yr and 5.6 cm at 40 yr for νa = 1017, 1018 and 1019 Pa s. For the case of νa = 1020, instead, the maximum geoid anomaly decreases at later times due to the fact that the positive contribution due to stress relaxation by viscous flow, U P(t), is smaller than the reduction, linearly in time, of the elastic response to the stick- and back-slips, U E(t). The short wavelengths of the elastic response remain visible at all times and make the maximum geoid anomaly always greater than that from the periodic contribution due to stress relaxation by viscous flow. This pattern is not altered even when the local minimum close to the trench between the two maxima appears at t = 2 and 20 yr for νa = 1017 and 1018 Pa s (Fig. 3). Furthermore, we note that the periodic geoid anomaly for νa = 1017 Pa s becomes smaller than that for νa = 1018, 1019 and 1020 Pa s between 2 to 4 yr, 4 to 10 yr and 10 to 20 yr. Particularly, at t = 200 yr which is one half the interseismic period, T/2, the amplitude of the geoid anomaly for νa = 1017 Pa s does not exceed 0.7 cm along the whole transect. At this time, instead, the geoid anomalies for νa = 1018, 1019 and 1020 Pa s reach maximum values of 1.6, 3.2 and 2.4 cm close to the trench. This results from the faster decay of the geoid anomaly due to stress relaxation for the lowest asthenospheric viscosity of 1017 Pa s. At later times, the geoid anomalies for all the asthenospheric viscosities decay to zero everywhere, as expected on the basis of the results obtained in Section 4 for the heuristic case of a discrete relaxation spectrum. In fact, the contribution from each relaxation mode reaches its maximum amplitude before one half the interseismic period, T/2 = 200 yr. In order to provide a clearer description of our results, we show in Fig. 5 the periodic geoid anomaly as function of time at two representative points on the side of the overriding plates, at x = −150 and −50 km. The point at x = −150 km is about where the maximum geoid anomaly during the early stage of the seismic cycle is located. The point at x = −50 km, instead, is about where the local minimum between the two maxima for the cases of νa = 1017 and 1018 Pa s is located. As expected, the periodic geoid anomaly due to only stress relaxation resembles that of the functions pj controlling the time evolution of the contributions of the relaxation modes, within the framework of the heuristic Earth model with a discrete relaxation spectrum. Differences, like the change in sign in Fig. 5(b) for νa = 1017 Pa s and in Fig. 5(d) for νa = 1017 and 1018 Pa s and the different concavities after the early phase of the seismic cycle for νa = 1017 and 1018 Pa s, are the result of the superimposition of relaxation modes with strengths of opposite signs. Figure 5. View largeDownload slide (a,c) Periodic geoid anomaly during one seismic cycle due to both the instantaneous elastic response and stress relaxation by viscous flow U (E(t) + P(t)); (b,d) as in (a,c) but limited to relaxation by viscous flow U P(t), at the distances from the trench x = −150 (a,b) and −50 km for (c,d). Solid, dashed, dashed-dotted and dotted lines refer to the case of asthenospheric viscosities of νa = 1017, 1018, 1019 and 1020 Pa s. The thick solid line in panels (a) and (c) indicates the instantaneous elastic response, U E(t). Figure 5. View largeDownload slide (a,c) Periodic geoid anomaly during one seismic cycle due to both the instantaneous elastic response and stress relaxation by viscous flow U (E(t) + P(t)); (b,d) as in (a,c) but limited to relaxation by viscous flow U P(t), at the distances from the trench x = −150 (a,b) and −50 km for (c,d). Solid, dashed, dashed-dotted and dotted lines refer to the case of asthenospheric viscosities of νa = 1017, 1018, 1019 and 1020 Pa s. The thick solid line in panels (a) and (c) indicates the instantaneous elastic response, U E(t). At the end, for the sake of completeness, we show in Fig. 6 also the static geoid anomaly from the steady-state response after the infinite series of seismic cycles, U S∞, on which the periodic change shown in Fig. 4 is superimposed. As anticipated at the end of Section 3, this contribution is smaller, with amplitudes of at most 3.2 cm, than the static geoid anomaly due to the long term average steady state motion at subduction zones and to the mantle density heterogeneities and the induced mantle flow, which is of the order of tens of meters (Matsu’ura & Sato 1989; Sato & Matsu’ura 1993; Zhong & Gurnis 1992; Marotta et al. 2006). Figure 6. View largeDownload slide The static contribution due to stress relaxation by viscous flow after an infinite number of seismic cycles, U S∞, for asthenospheric viscosities of νa = 1017, 1018, 1019 and 1020 Pa s (solid, dashed, dashed-dotted, and dotted lines). Figure 6. View largeDownload slide The static contribution due to stress relaxation by viscous flow after an infinite number of seismic cycles, U S∞, for asthenospheric viscosities of νa = 1017, 1018, 1019 and 1020 Pa s (solid, dashed, dashed-dotted, and dotted lines). 6 CONCLUSION We have extended the approach proposed by Savage (1983) from gravitating and flat to self-gravitating and spherically symmetric Earth models in order to quantify the geoid anomaly caused by the strain accumulation and release at subduction zones. Different from Savage (1983) that considered periodic topography changes during the seismic cycles, we have focused on the geoid anomaly and obtained analytical expressions for the periodic geoid anomaly after an infinite number of seismic cycles, as well as for the a smaller static contribution from the steady-state response over which this periodic change is superimposed. These analytical expressions do not require to calculate the viscoelastic response to each earthquake in the time domain, as done instead by Savage (1983) on the basis of the solution provided by Thatcher & Rundle (1981). Instead, they are based on the Correspondence Principle between the viscoelastic problem in the time domain and an equivalent elastic problem in the Laplace domain, and the analysis of the singularities of the geoid anomaly in the Laplace domain discussed in Appendix B. Further, within the heuristic assumption of a Maxwell Earth model with a discrete relaxation spectrum, we have understood how the contribution to the periodic geoid anomaly from relaxation modes with characteristic relaxation times comparable or greater than the interseismic period is damped, until it is completely filtered out when the relaxation times are greater than the interseismic period. We have applied this methodology to the case of an infinite sequence of megathrust earthquakes with moment magnitude Mw = 9 and an interseismic period of T = 400 yr, and computed the periodic geoid anomaly during the seismic cycle. The elastic response at the earthquake time yields a positive geoid anomaly with a maximum value of 3.8 cm close to the trench and on the side of the overriding plate and a negative geoid anomaly further away from the trench with a minimum value of −1.3 cm (Fig. 4 at t = 0 yr). At later times, this instantaneous elastic response decreases to zero linearly in time due to the elastic contribution associated with the locking of the fault during the interseismic phase. During the interseismic phase, the periodic geoid anomaly is affected also by stress relaxation by viscous flow. The latter provides a peculiar pattern of global upwarping of the geoid that is broader in space compared to the elastic contribution and has its maximum close to the trench and on the side of the overriding plate (Fig. 3). This positive geoid anomaly grows in the early stage up to values that can be two times greater than the maximum amplitude of the elastic response to the earthquakes for low asthenospheric viscosities. Particularly, the maximum geoid anomalies due to stress relaxation attain maximum values of 5.2 at 1 and 10 yr after the earthquake for asthenospheric viscosities of 1017 and 1018 Pa s. In contrast, for the case of an asthenospheric viscosities of 1019 and 1020 Pa s, the maximum positive geoid anomalies attain maximum values of 4.2 and 0.9 cm at 70 and 150 yr after the earthquake. This is consistent with the finding that the interaction between stress relaxation by viscous flow and the periodic dip-slip history describing strain accumulation and release damps the contribution from those relaxation modes which have relaxation times comparable or greater than the interseismic period (Fig. 2). Once the geoid anomaly due to stress relaxation has attained its maximum, it starts to decrease to zero. In this second phase of the seismic cycle, we note that the lowermost asthenospheric viscosity of 1017 Pa s yields the fastest collapse of the geoid maximum which is splitted into two reduced maxima intermingled by a local minimum close to the trench, finally developing into a negative geoid anomaly starting from t = 10 yr after the earthquake. These results constitutes a first step towards the interpretation of time-dependent geoid anomalies in the surroundings of subduction zones and can be used as a theoretical framework for explaining data time series from the GRACE and GOCE satellite missions even at those locations where no megathurst earthquakes have occurred during the time in which these two missions have been operative, from 2002 to present and from 2009 to 2013. Particularly, an improved time series analysis by Meyer et al. (2016), based on the optimization of the GRACE data, leads to the new release AIUB (Astronomical Institute of the University of Bern) — RL02 and allows to estimate the standard deviation (formal errors) for the time-dependent geoid evaluation up to the harmonic degree ℓ = 60. This standard deviation, evaluated from the difference between the time-dependent geoid and a static reference geoid, which is concordant with our approach where the time-dependent geoid must be considered as superimposed to the static geoid due to steady-state subduction, is about 1 mm (fig. 2 of Meyer et al. 2016) and indicates that our co-, post- and interseismic signals are detectable from GRACE, since in the interseismic phase the geoid anomaly grows to several centimetres on timescales ranging from one year to a few decades, depending on the asthenospheric viscosities. Although our results take into account all the wavelengths of the geoid anomaly, in view of its large spatial scale, the filtering in the bandwidth up to the harmonic degree ℓ = 60 of GRACE does not affect this conclusion. Indeed, this filtering reduces the amplitude of the geoid anomalies due to the elastic response, U E, and to stress relaxation by viscous flow, U P, by about a factor of 4 and 2, respectively, with the largest reduction affecting the elastic response because of its shorter wavelengths. The comparison of our modelling with space gravity data, in addition to proper wavelength filtering, requires specific assumptions on the locked portion of the fault during the interseismic phase, which does not necessarily coincide with that adopted in this work (extending form the Earth surface down to the outermost 50 km), and the knowledge of past seismicity back in time as far as the post-seismic contribution significantly affects the data time series. The sensitivity analysis to different locking of the fault and to the realistic earthquake sequences will be the subject of future works. Acknowledgements This research was partially supported by the CAS/CAFEA international partnership program for creative research teams (No. KZZDEW-TZ-19). DY acknowledges support from Geochemistry and CISE from National Science Foundation. REFERENCES Broerse T., Riva R., Vermeersen B., 2011. Ocean contribution to seismic gravity changes: the sea level equation for seismic perturbations revisited, Geophys. J. Int. , 199, 1094– 1109. Google Scholar CrossRef Search ADS   Broerse T., Riva R., Simons W., Govers R., Vermeersen B., 2015. Postseismic grace and GPS observations indicate a rheology contrast above and below the Sumatra slab, J. geophys. Res. , 120, 5343– 5361. Google Scholar CrossRef Search ADS   Cambiotti G., Sabadini R., 2010. The compressional and compositional stratifications in Maxwell Earth models: the gravitational overturning and the long-period tangential flux, Geophys. J. Int. , 180, 475– 500. Google Scholar CrossRef Search ADS   Cambiotti G., Sabadini R., 2013. Gravitational seismology retrieving centroid moment tensor solution of the 2011 Tohoku earthquake, J. geophys. Res. , 118, 183– 194. Google Scholar CrossRef Search ADS   Cambiotti G., Sabadini R., 2015. On the seismic perturbation due to a fault system: its evaluation beyond the epicentral reference frame, Geophys. J. Int. , 203, 943– 959. Google Scholar CrossRef Search ADS   Cambiotti G., Bordoni A., Sabadini R., Colli L., 2011a. Grace gravity data help constraining seismic models of the 2004 Sumatran earthquake, J. geophys. Res. , 116, B10403. Google Scholar CrossRef Search ADS   Cambiotti G., Ricard Y., Sabadini R., 2011b. New insights into mantle convection true polar wander and rotational bulge readjustment, Earth planet. Sci. Lett. , 310, 538– 543. Google Scholar CrossRef Search ADS   Cambiotti G. Y. R., Sabadini R., 2011c. Ice age True PolarWander in a compressible and non-hydrostatic Earth, Geophys. J. Int. , 183, 1248– 1264. Google Scholar CrossRef Search ADS   Cambiotti G., Klemann V., Sabadini R., 2013. Compressible viscoelastodynamics of a spherical body at long timescales and its isostatic equilibrium, Geophys. J. Int. , 193, 1071– 1082. Google Scholar CrossRef Search ADS   Chen J., Wilson C., Tapley B., Grand S., 2007. Grace detects coseismic and postseismic deformation from the sumatra-andaman earthquake, Geophys. Res. Lett. , 34, L13302. De Linage C., Rivera L., Hinderer J., Boy J.-P., Rogister Y., Lambotte S., Biancale R., 2009. Separation of coseismic and postseismic gravity changes for the 2004 Sumatra–Andaman earthquake from 4.6 yr of GRACE observations and modelling of the coseismic change by normal-modes summation, Geophys. J. Int. , 176, 695– 714. Google Scholar CrossRef Search ADS   Dziewonski A., Anderson D., 1981. Preliminary Reference Earth Model, Phys. Earth planet. Inter. , 25, 297– 356. Google Scholar CrossRef Search ADS   Fuchs M., Hooper A., Broerse T., Bouman J., 2016. Distributed fault slip model for the 2011 Tohoku-Oki earthquake from GNSS and GRACE/GOCE satellite gravimetry, J. geophys. Res. , 121, 1114– 1130. Google Scholar CrossRef Search ADS   Gross R., Chao B., 2001. The gravitational signature of earthquakes, in Gravity, Geoid and Geodynamics 2000 , 123, pp. 205– 210, ed. Sideris M.G., International Association of Geodesy Symposia. Google Scholar CrossRef Search ADS   Han S., Sauber J., Luthcke S., Ji C., Pollitz F., 2008. Implications of postseismic gravity change following the great 2004 Sumatra–Andaman earthquake from the regional harmonic analysis of grace intersatellite tracking data, J. geophys. Res. , 113, B11413. Google Scholar CrossRef Search ADS   Marotta A., Spelta E., Rizzetto C., 2006. Gravity signature of crustal subduction inferred from numerical modelling, Geophys. J. Int. , 166, 923– 938. Google Scholar CrossRef Search ADS   Matsuo K., Heki K., 2011. Coseismic gravity changes of the 2011 Tohoku-Oki earthquake from satellite gravimetry, Geophys. Res. Lett. , 38, L00G12. Google Scholar CrossRef Search ADS   Matsu’ura M., Sato T., 1989. A dislocation model for the earthquake cycle at convergent plate boundaries, Geophys. J. Int. , 96, 23– 32. Google Scholar CrossRef Search ADS   Meyer U., Jäggi A., Jean Y., Beutler G., 2016. Aiub-rl02: an improved time-series of monthly gravity fields from grace data, Geophys. J. Int. , 205, 1196– 1207. Google Scholar CrossRef Search ADS   Plat H.-P., Jüttner H.-U., 1995. Rayleigh-Taylor instabilities of a self-gravitating Earth, J. Geodyn. , 20, 267– 288. Google Scholar CrossRef Search ADS   Pollitz F., Bürgmann R., Banerjee P., 2006. Post-seismic relaxation following the great 2004 Sumatra–Andaman earthquake on a compressible self-gravitating earth, Geophys. J. Int. , 167, 397– 420. Google Scholar CrossRef Search ADS   Ricard Y., Richards M., Lithgow-Bertelloni C., Lestunff Y., 1993. A geodynamic model of mantle density heterogeneity, J. geophys. Res. , 98, 21 895– 21 909. Google Scholar CrossRef Search ADS   Sato T., Matsu’ura M., 1993. A kinematic model for evolution of island arc-trench systems, Geophys. J. Int. , 114, 512– 530. Google Scholar CrossRef Search ADS   Savage J.C., 1983. A dislocation mode of strain accumulation and release at a subduction zone, J. geophys. Res. , 88, 4984– 4996. Google Scholar CrossRef Search ADS   Savage J.C., Burford R.O., 1973. Geodetic determination of relative plate motion in Central California, J. geophys. Res. , 78, 832– 845. Google Scholar CrossRef Search ADS   Segall P., 2010. Earthquake and Volcano Deformation , Princeton University Press. Google Scholar CrossRef Search ADS   Song T., Simons M., 2003. Large Trench-Parallel Gravity Variations Predict Seismogenic Behavior in Subduction Zones, Science , 301, 630– 633. Google Scholar CrossRef Search ADS PubMed  Sun W., Okubo S., 2004. Co-seismic deformations detectable by satellite gravity missions, J. geophys. Res. , 109, B04405. Tanaka Y., Okuno J., Okubo S., 2006. A new method for the computation of global viscoelastic post-seismic deformation in a realistic earth model (I)–vertical displacement and gravity variation, Geophys. J. Int. , 164, 273– 289. Google Scholar CrossRef Search ADS   Tanaka Y., Hasegawa T., Tsuruoka H., Klemann V., Martinec Z., 2015. Spectral-finite element approach to post-seismic relaxation in a spherical compressible Earth: application to gravity changes due to the 2004 Sumatra–Andaman earthquake, Geophys. J. Int. , 200, 299– 321. Google Scholar CrossRef Search ADS   Thatcher W., Rundle J., 1979. A model for the earthquake cycle in underthrust zones, J. geophys. Res. , 84, 5540– 5556. Google Scholar CrossRef Search ADS   Thatcher W., Rundle J., 1981. A viscoelastic model for periodically recurring earthquakes in subduction zones (abstract), EOS, Trans. Am. geophys. Un. , 62, 1031. Tosi N., Sabadini R., Marotta A. M., Vermeersen L., 2005. Simultaneous inversion for the Earth’s mantle viscosity and ice mass imbalance in Antarctica and Greenland, J. geophys. Res. , 110, B07402. Google Scholar CrossRef Search ADS   Wang L., Shum C.F.J.S., Tapley B., Dai C., 2012. Coseismic and postseismic deformation of the 2011 Tohoku-Oki earthquake constrained by GRACE gravimetry, Geophys. Res. Lett. , 39, L07301. Wolf D., 1991. Visco-elastodynamics of a stratified compressible planet: incremental field equations and short- and long-time asymptotes, Geophys. J. Int. , 104, 401– 417. Google Scholar CrossRef Search ADS   Wolf D., Kaufmann G., 2000. Effects due to compressional and compositional density stratification on load-induced Maxwell-viscoelastic perturbations, Geophys. J. Int. , 140, 51– 62. Google Scholar CrossRef Search ADS   Zhong S., Gurnis M., 1992. Viscous flow model of a subduction zone with a faulted lithosphere: Long and short wavelength topography, gravity, and geoid, Geophys. Res. Lett. , 19, 1891– 1894. Google Scholar CrossRef Search ADS   APPENDIX A: ELASTIC AND VISCOELASTIC GREEN’s FUNCTIONS Within the framework of self-gravitating spherically symmetric elastic Earth models and adapting the formalism presented in Cambiotti & Sabadini (2015) for the case of a time-dependent uniform slip over the fault surface (i.e. the slip varies through time but it does not depend on the position over the fault), the geoid anomaly at a point of the Earth surface can be expressed as follows:   \begin{eqnarray} G(\theta ,\varphi ,t) = \sum _{j=1}^{2} \delta u_j(t)\, K_E^j(\theta ,\varphi ) \end{eqnarray} (A1)where δuj are the along strike (j = 1) and along dip (j = 2) components of the uniform slip, and $$K^j_E(\theta ,\varphi )$$ are the elastic Green functions for the geoid anomaly at colatitude θ and longitude φ given by   \begin{eqnarray} K^j_E(\theta ,\varphi )=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell } K_{E\ell m}^j\, Y_{\ell m}(\theta ,\varphi ). \end{eqnarray} (A2)Here, Yℓm are the normalized real spherical harmonics of degree ℓ and order m, and $$K_{E\ell m}^j$$ are the spherical harmonic coefficients of $$K^j_E(\theta ,\varphi )$$. The latter can be obtained as follows:   \begin{eqnarray} K_{E\ell m}^j = a^{-2}\, \sum _{p=0}^{2}\int _{\mathcal {S}} k_{\ell ,p}(a,r^{\prime })\, \Delta _{\ell m,pj}^D(\mathbf {r}^{\prime })\, \mathrm{d}S^{\prime } \end{eqnarray} (A3)where a is the Earth radius, $$\mathcal {S}$$ is the fault surface, dS΄ is the infinitesimal area, r΄ is the vector position identifying a point of the fault surface, kℓ, p(a, r΄) is the (elastic) gravitational seismic Love number governing the polar (p = 0), bipolar (p = 1) and quadrupolar (p = 2) patterns of the response at the Earth surface to a point-like seismic force at the radial distance from the Earth centre r΄ = ‖r΄‖, and Δℓm, pj is the following factor:   \begin{eqnarray} \Delta _{\ell m,pj}(\mathbf {r}^{\prime })= \frac{D_j^p(\mathbf {r}^{\prime })\, F_{\ell m}^p(\theta ^{\prime },\varphi ^{\prime })+D_j^{-p}(\mathbf {r}^{\prime })\, F_{\ell m}^{-p}(\theta ^{\prime },\varphi ^{\prime })}{1+\delta _{p0}} \end{eqnarray} (A4)which accounts for the geometry (dip and strike angles) and location (colatitude θ΄ and longitude φ΄) of the infinitesimal surface element, with $$F_{\ell m}^p$$ and $$D_{j}^p$$ given by eqs (31) and (33) of Cambiotti & Sabadini (2015). In order to perform the integration over depth z΄ (or the radial distance from the Earth surface r΄ = a − z΄) we interpolate the seismic Love numbers sampled every 1 km, but for the shallowest 3 km where the sampling is gradually refined up to 0.1 km. In the main text we will not explicit the dependence on the colatitude and longitude of the point at which the geoid anomaly is computed and, as far as we will focus only on dip-slip faults, we will omit the indexing for distinguishing between strike-slip and dip-slip and we simply write   \begin{eqnarray} G(t) = K_E\, \delta u(t) \end{eqnarray} (A5)where it is understood that δu and KE correspond to δu2 and $$K^2_E$$ in the formulation used in this appendix. According to the correspondence principle (Tanaka et al. 2006), these results can be generalized to Maxwell Earth models as follows:   \begin{eqnarray} G(t) = (K\star \delta u)(t) \end{eqnarray} (A6)denoting the time convolution between the slip time history δu and the viscoelastic Green function K. The latter is defined as the inverse Laplace transform of the equivalent elastic Green functions $$\tilde{K}$$ in the Laplace domain   \begin{eqnarray} K(t)=\mathcal {L}^{-1}[\tilde{K}(s)] \end{eqnarray} (A7)where $$\tilde{K}(s)$$ is obtained solving the equivalent elastic problem where the shear modulus μ is replaced by the equivalent shear modulus $$\hat{\mu }(s),$$  \begin{eqnarray} \hat{\mu }(s) = \frac{\mu s}{s+1/\tau }, \end{eqnarray} (A8)providing the dependence on the Laplace variable s. The lambda parameter λ, instead, is replaced by   \begin{eqnarray} \hat{\lambda }(s) = \kappa - \frac{2}{3}\, \hat{\mu }(s). \end{eqnarray} (A9) In order to distinguish between the instantaneous elastic response and the contributions caused by stress relaxation by viscous flow, the viscoelastic Green function is usually decomposed as follows:   \begin{eqnarray} K(t) = K_E\, \delta (t)+K_V(t) \end{eqnarray} (A10)with   \begin{eqnarray} K_V(t)=\mathcal {L}^{-1}[\tilde{K}_V(s)] \end{eqnarray} (A11a)  \begin{eqnarray} \tilde{K}_V(s) = \tilde{K}(s)-K_E. \end{eqnarray} (A11b) Then, substituting eq. (A10) into eq. (A6), we obtain   \begin{eqnarray} G(t) = G_E(t)+G_V(t), \end{eqnarray} (A12)where GE and Gv are the instantaneous elastic response and the response due to only stress relaxation by viscous flow,   \begin{eqnarray} G_E(t) =\delta u(t)\, K_E \end{eqnarray} (A13a)  \begin{eqnarray} G_V(t) =\mathcal {L}^{-1}[\tilde{G}_V(s)] \end{eqnarray} (A13b)with   \begin{eqnarray} \tilde{G}_V(s)=\tilde{K}_V(s)\, \delta \tilde{u}(s). \end{eqnarray} (A14)Here, we have made use of the fact that the time convolution between two functions corresponds to the product of their Laplace transforms in the Laplace domain. In the summation over the harmonic degrees in eq. (A2) we have considered the harmonic degrees up to 5 × 105 and 2 × 103 in order to obtain the elastic, GE, and viscous, Gv, Green functions, respectively. We note that the viscous Green function requires a much lower truncation than the elastic one because it reflects the stress relaxation by viscous flow that occurs, within the assumption of the present work, at depths greater than 50 km. APPENDIX B: GEOID ANOMALY In order to estimate the geoid anomaly due to only stress relaxation by viscous flow, we first consider the Laplace transform of the dip-slip time history that, by making use of the expression for geometric series, we rearrange as follows:   \begin{eqnarray} \delta \tilde{u}(s) = U\, \left(-\frac{1}{s^2\, T}+\sum _{n=1}^{\infty } \frac{e^{-n\, s\, T}}{s}\right)=U\, \frac{\tilde{g}(s)}{s} \end{eqnarray} (B1)with   \begin{eqnarray} \tilde{g}(s) = -\frac{1}{1-e^{s\, T}}-\frac{1}{s\, T} \end{eqnarray} (B2)Then, substituting eq. (B1) into eq. (A14), we obtain   \begin{eqnarray} \tilde{G}_V(s) = U\, \frac{\tilde{K}_V(s)\, \tilde{g}(s)}{s}. \end{eqnarray} (B3) In order to perform the inverse Laplace transform of $$\tilde{G}_V(s)$$, we must first investigate its singularities in the Laplace domain. They are depicted in Fig. B1. The function $$\tilde{K}_V(s)$$ has singularities on the negative half of the real axis of the Laplace domain and, particularly, they are confined to the interval $$(-1/\tau _{\rm {min}})\subset \mathbb {R}^-$$. These singularities depend on the specific viscoelastic stratification of the Earth model and can be isolated poles, clusters of poles and continuous not-Lipschitzian intervals (Tanaka et al. 2006; Cambiotti et al. 2011c). The function $$\tilde{g}(s)$$, instead, has singularities at the isolated poles s = i ωn (n = ±1, …, ±∞) with   \begin{eqnarray} \omega _n = \frac{2\, \pi \, n}{T} \end{eqnarray} (B4)and its residues are   \begin{eqnarray} \lim _{s\rightarrow \infty } \tilde{g}(s)\, (s-i\, \omega _n) = \frac{1}{T}. \end{eqnarray} (B5)We note that s = 0 is not a singularity for $$\tilde{g}(s)$$ and, in fact, its limit for s → 0 exists   \begin{eqnarray} \lim _{s\rightarrow 0} \tilde{g}(s) = -\frac{1}{2}. \end{eqnarray} (B6) Figure B1. View largeDownload slide Schematic depicting the singularities of the function $$\tilde{K}_V(s)\, \tilde{g}(s)$$ in the Laplace domain. The black circles indicate the poles along the imaginary axis due to $$\tilde{g}(s)$$, which are at s = i ωn (n = ±1, …, ±∞). The white circles and rectangles indicate the poles and the not-Lipschitzian intervals due to $$\tilde{K}_V(s)$$. The latter lay on the negative half of the real axis of the Laplace domain and, particularly, within the interval $$(-1/\tau _{\rm {min}},0)\subset \mathbb {R}^-$$ identified by the thick solid line. The dashed rectangle is the closed contour Γ that includes the origin of the Laplace domain and all the singularities due to $$\tilde{K}_V(s)$$. Figure B1. View largeDownload slide Schematic depicting the singularities of the function $$\tilde{K}_V(s)\, \tilde{g}(s)$$ in the Laplace domain. The black circles indicate the poles along the imaginary axis due to $$\tilde{g}(s)$$, which are at s = i ωn (n = ±1, …, ±∞). The white circles and rectangles indicate the poles and the not-Lipschitzian intervals due to $$\tilde{K}_V(s)$$. The latter lay on the negative half of the real axis of the Laplace domain and, particularly, within the interval $$(-1/\tau _{\rm {min}},0)\subset \mathbb {R}^-$$ identified by the thick solid line. The dashed rectangle is the closed contour Γ that includes the origin of the Laplace domain and all the singularities due to $$\tilde{K}_V(s)$$. We can take advantage of the fact that we have analytical expressions for the poles and the residues of $$\tilde{g}(s)$$ and that $$\tilde{K}_V(s)$$ is analytical at these poles decomposing the product $$\tilde{K}_V(s)\, \tilde{g}(s)$$ as follows:   \begin{eqnarray} \tilde{K}_V(s)\, \tilde{g}(s) = \tilde{p}(s)+\tilde{c}(s), \end{eqnarray} (B7)where $$\tilde{p}(s)$$ accounts for all the poles at s = i ωn:   \begin{eqnarray} \tilde{p}(s) = \frac{1}{T}\, \sum _{n=1}^\infty \left(\frac{\tilde{K}_V(i\, \omega _n)}{s-i\, \omega _n}+ \frac{\tilde{K}_V^*(i\, \omega _{n})}{s+i\, \omega _n}\right). \end{eqnarray} (B8)Here * stands for the complex conjugate and we have made use of the fact that the Laplace transforms of real functions has the following property: $$K_V(s^*) = K_V^*(s)$$. The function $$\tilde{c}(s)$$, instead, is simply the difference between the product $$\tilde{K}^V(s)\, \tilde{g}(s)$$ and $$\tilde{p}(s),$$  \begin{eqnarray} \tilde{c}(s) = \tilde{K}_V(s)\, \tilde{g}(s)-\tilde{p}(s). \end{eqnarray} (B9)This function has the same singularities of $$\tilde{K}_V(s)$$, but it has not those from $$\tilde{g}(s)$$ because they are removed by subtracting $$\tilde{P}(s)$$. Indeed, the residues of $$\tilde{C}(s)$$ at s = i ωn are zero. Within this decomposition, the geoid anomaly due to stress relaxation by viscous flow, eq. (A13b), that is the inverse Laplace transform of eq. (B3), can be written as follows:   \begin{eqnarray} G_V(t) = U \left(P(t)+C(t)\right)\!, \end{eqnarray} (B10)where P(t) is the inverse Laplace transform of $$\tilde{p}(s)/s$$, which reads   \begin{eqnarray} P(t) \!=\! \mathcal {L}^{-1}\big [\tilde{p}(s)/s\big ] \!=\! \sum _{n=1}^\infty \big [I_n\, (1\!-\!\cos (\omega _n\, t))\!+\!R_n\, \sin (\omega _n\, t)\big ] \end{eqnarray} (B11)with   \begin{eqnarray} R_n = \frac{\Re [K_V(i\, \omega _n)]}{\pi \, n} \end{eqnarray} (B12a)  \begin{eqnarray} I_n = \frac{\Im [K_V(i\, \omega _n)]}{\pi \, n} \end{eqnarray} (B12b)and C(t) is the inverse Laplace transform of $$\tilde{c}(s)/s$$, which can be obtained by integration in the Laplace domain over a closed contour Γ, containing all its singularities,   \begin{eqnarray} C(t) = \mathcal {L}^{-1}[\tilde{c}(s)/s] = \frac{1}{2\, \pi \, i}\, \int _{\Gamma } \frac{\tilde{c}(s)}{s}\, e^{s\, t}\, \mathrm{d}s. \end{eqnarray} (B13)We note that, as depicted in Fig. B1, the closed contour Γ must contain the origin s = 0 and the interval $$(-1/\tau _{\rm {min}},0)\subset \mathbb {R}^-$$ of the negative half of the real axis of the Laplace domain. As implied by the sine and cosine functions in eq. (B11), P(t) is a periodic function with period T and yields zero at t = n T. On the contrary, due to the singularities of $$\tilde{c}(s)$$ having negative real part and zero imaginary part, C(t) has an asymptotic limit for t → ∞, say C∞, that we can estimate by taking the limit of $$\tilde{c}(s)$$ for s → ∞   \begin{eqnarray} C_\infty &=& \lim _{t\rightarrow \infty } C(t) = \lim _{s\rightarrow 0} \tilde{c}(s) = -\frac{\tilde{K}_V(0)}{2} - \tilde{p}(0) \nonumber \\ &=& -\frac{\tilde{K}_V(0)}{2}+\sum _{n=1}^{\infty } I_n \end{eqnarray} (B14)by making use of eqs (B6), (B8) and (B9). APPENDIX C: DISCRETE RELAXATION SPECTRUM The expression obtained in Appendix A greatly simplifies when we consider the case in which the relaxation spectrum of the Maxwell Earth model is discrete. In this case, the viscoelastic Green function takes the following form:   \begin{eqnarray} K(t) = K_E\, \delta (t)+ H(t)\, \sum _{j}k_j\, e^{-t/T_j}\!, \end{eqnarray} (C1)where δ(t) is the Dirac delta, and kj and Tj are the residue and the characteristic relaxation time of the jth relaxation mode. This form, in the Laplace domain, becomes   \begin{eqnarray} \tilde{K}(s) = K_E+\sum _j\frac{k_j}{s-s_j} \end{eqnarray} (C2)with sj = −1/Tj being the pole of the jth relaxation mode, confined on the interval $$[-1/\tau _{\rm {min}},0]\subset \mathbb {R}^-$$ of the negative half of the real axis of the Laplace domain. In this case and considering the dip-slip time history describing the infinite earthquake sequence and back-slip, making use of the residue theorem, the geoid anomaly becomes   \begin{eqnarray} {G(t)= U\, \sum _j\bigg [K_E\, f(t)+r_j\, \left(\frac{1-e^{-t/T_j}}{T/T_j}-\frac{t}{T}\right)} \nonumber \\ {\qquad\qquad\qquad+ \sum _{n=1}^{\infty } r_j\, \big (1-e^{-(t-n\, T)/T_j}\big )\, H(t-n\, T)\bigg ] } \end{eqnarray} (C3)where rj = −kj/sj = kj Tj is the so called strength of the jth relaxation mode, that is the contribution of the mode to the geoid anomaly due to a stick-slip earthquake in the fluid limit, for t → ∞. The first term in the square brackets of the RHS of eq. (C3) is the instantaneous elastic response, the second term is the response to the back-slip due to only stress relaxation and the third term is the summation of the post-seismic response to each earthquake. Considering a time during the nth seismic cycle, t + n T with t ∈ [0, T), and using the expression for finite geometric series for dealing with the post-seismic contribution of the first nth earthquakes, after some algebra, the geoid anomaly during the nth seismic cycle, eq. (10), becomes   \begin{eqnarray} {G_n(t) = U\, \bigg [-K_E\, \frac{t}{T} +\sum _j r_j \bigg (\frac{1-e^{-(t+n\, T)/T_j}}{T/T_j} -\frac{t}{T}} \nonumber \\ {\qquad\qquad\qquad- e^{(T-t)/T_j}\, \frac{1-e^{-n\, T/T_j}}{1-e^{-T/T_j}}\bigg )\bigg ].} \end{eqnarray} (C4)Then, reorganizing this result as in eq. (10), we obtain specific analytical expressions for P(t) and Sn(t), with E(t) given by eq. (11a). Particularly, we have   \begin{eqnarray} P(t) = \sum _jr_j\, p_j(t) \end{eqnarray} (C5a)  \begin{eqnarray} S_n(t) = \sum _j r_j \left(1-e^{-(n\, T+t)/T_j}\right)\left(\frac{T_j}{T}-\frac{1}{1-e^{-T/T_j}}\right) - K_E \end{eqnarray} (C5b)with   \begin{eqnarray} p_j(t) = \frac{1-e^{s_j\, t}}{1-e^{s_j\, T}}-\frac{t}{T}. \end{eqnarray} (C6)We note that   \begin{eqnarray} S_\infty = \sum _j r_j\, \left(\frac{T_j}{T}-\frac{1}{1-e^{-T/T_j}}\right)-K_E. \end{eqnarray} (C7) It can be verified that these expressions are consistent with eqs (B11) and (B14), considering that, in view of eq. (C2), eq. (B12) becomes   \begin{eqnarray} R_n = \frac{1}{\pi \, n}\sum _j\frac{r_j}{1+(T_j\, \omega _n)^2} \end{eqnarray} (C8a)  \begin{eqnarray} I_n = -\sum _j\frac{2\, r_j\, T_j/T}{1+(T_j\, \omega _n^2)} \end{eqnarray} (C8b)and we have   \begin{eqnarray} \tilde{K}_V(0) = \sum _j r_j. \end{eqnarray} (C9) © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Time-dependent geoid anomalies at subduction zones due to the seismic cycle

, Volume 212 (1) – Jan 1, 2018
12 pages

/lp/ou_press/time-dependent-geoid-anomalies-at-subduction-zones-due-to-the-seismic-Hq1WfUYcaU
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx421
Publisher site
See Article on Publisher Site

### Abstract

Abstract We model the geoid anomalies excited during a megathrust earthquake cycle at subduction zones, including the interseismic phase and the contribution from the infinite series of previous earthquakes, within the frame of self-gravitating, spherically symmetric, compressible, viscoelastic Earth models. The fault cuts the whole 50 km lithosphere, dips 20°, and the slip amplitude, together with the length of the fault, are chosen in order to simulate an Mw = 9.0 earthquake, while the viscosity of the 170 km thick asthenosphere ranges from 1017 to 1020 Pa s. On the basis of a new analysis from the Correspondence Principle, we show that the geoid anomaly is characterized by a periodic anomaly due to the elastic and viscous contribution from past earthquakes and to the back-slip of the interseismic phase, and by a smaller static contribution from the steady-state response to the previous infinite earthquake cycles. For asthenospheric viscosities from 1017–1018 to 1019–1020 Pa s, the characteristic relaxation times of the Earth model change from shorter to longer timescales compared to the 400 yr earthquake recurrence time, which dampen the geoid anomaly for the higher asthenospheric viscosities, since the slower relaxation cannot contribute its whole strength within the interseismic cycle. The geoid anomaly pattern is characterized by a global, time-dependent positive upwarping of the geoid topography, involving the whole hanging wall and partially the footwall compared to the sharper elastic contribution, attaining, for a moment magnitude Mw = 9.0, amplitudes as high as 6.6 cm for the lowermost asthenospheric viscosities during the viscoelastic response compared to the elastic maximum of 3.8 cm. The geoid anomaly vanishes due to the back-slip of the interseismic phase, leading to its disappearance at the end of the cycle before the next earthquake. Our results are of importance for understanding the post-seismic and interseismic geoid patterns at subduction zones. Seismic cycle, Subduction zone processes, Dynamics: gravity and tectonics, Dynamics: seismoctectonics 1 INTRODUCTION The 2004 Sumatra, 2010 Maule and 2011 Tohoku megathrust earthquakes have caused the flourishing of studies on the gravity signatures caused by the earthquake-induced mass redistribution, nowadays detectable by GRACE and GOCE satellite missions. The co-seismic gravity anomalies caused by earthquakes with moment magnitude greater than 8.5 have been examined by means of long-wavelength gravity data from GRACE (Gross & Chao 2001; Sun & Okubo 2004; Matsuo & Heki 2011) and explained in terms of mass redistribution, volume variations and topography perturbations by De Linage et al. (2009), Cambiotti et al. (2011a) and Broerse et al. (2011). For the case of the 2011 Tohoku earthquake, GRACE data have been used for estimating its focal mechanism and hypocentre (Wang et al. 2012; Cambiotti & Sabadini 2013), and, jointly with GPS and GOCE data, to improve the constraint of its slip distribution over fault surface (Fuchs et al. 2016). In addition to the co-seismic signature, the post-seismic response associated to stress relaxation by viscous flow of the 2004 Sumatra earthquake has been investigated in a series of work (Pollitz et al. 2006; Chen et al. 2007; Han et al. 2008; Broerse et al. 2015; Tanaka et al. 2015) with the aim of constraining the viscosity stratification of the lithosphere and of the asthenosphere. Due to the monthly resolution of the GRACE satellite mission, estimates of the co-seismic jump at the earthquake time in the data time series are usually performed jointly with the estimate of the post-seismic response after the earthquake. Nevertheless, the interpretation of the latter requires a full understanding of all the contributions from the seismic cycle. Indeed, the time-dependent signature after the earthquake, in addition to the post-seismic response to the earthquake, is affected also by stress relaxation by viscous flow from past seismicity, as well as by strain accumulation due to the locking of the fault during the interseismic phase. As it concerns the topography changes at subduction zones, this issue has been first addressed by Savage (1983), accounting for both strain accumulation and release within the frame of viscoelastic, flat Earth models (Thatcher & Rundle 1979). According to Savage & Burford (1973) and Savage (1983) and as shown by Matsu’ura & Sato (1989) and Sato & Matsu’ura (1993) for old subduction zones, the long term average motion between the overriding and subducting plates is responsible for static vertical displacements, as for the peculiar bipolar gravity, geoid and topography patterns at subduction zones (Zhong & Gurnis 1992; Song & Simons 2003; Marotta et al. 2006), but it does not contribute to observable changes through time. In addition to the instantaneous elastic response at the earthquake time, changes in time during the interseismic phase are the results of the post-seismic response to the last earthquake, as well as to all the past seismicity and to the viscoelastic response to strain accumulation. Only for young subduction zones (i.e. for timescales smaller than the characteristic relaxation time of the high viscosity lithosphere) this relative motion yields a growth rate of the topography of island arc-trench systems. In this work, we will consider geoid anomalies at subduction zones due to both seismic and interseismic processes. Particularly, in order to simulate the geoid anomalies due to the locking of the fault during the interseismic phase, we follow the approach first proposed by Savage & Burford (1973) for strike-slip faults and, then, generalized to dip-slip faults in Savage (1983) (see also Segall 2010). According to this approach, the strain accumulation due to the locking of the fault is taken into account by imposing a constant in time slip rate on the fault system associated to the plate boundary that is opposite to the average slip rate necessary to accommodate the relative motion between plates. Compared to Savage (1983), we implement a compressible viscoelastic Earth model which is spherically symmetric and self-gravitating, thus allowing for a self-consistent treatment of the geoid anomalies, and builds on the elastic and density stratifications of the preliminary reference Earth model (PREM, Dziewonski & Anderson 1981). Furthermore, rather than considering separately the response to each earthquake and to the strain accumulation, we take advantage of the Correspondence Principle between the viscoelastic problem in the time domain and the equivalent elastic problem in the Laplace domain (Tanaka et al. 2006; Cambiotti et al. 2011c) in order to evaluate analytically the static geoid anomaly over which the periodic changes during the interseismic phase are superimposed. 2 SLIP TIME HISTORY As shown in Fig. 1, we approximate the fault surface at the subduction zone with a planar fault dipping with an angle α = 20° and we assume that the locked portion during the interseismic phase extends through the whole 50 km thick outermost elastic layer of the Earth model. Figure 1. View largeDownload slide Schematicdepicting the 50 km thick outermost elastic layer, the 170 km thick asthenosphere and upper and lower mantle below. The locked portion of the fault, indicated by the dashed line, dips with angle α = 20° and extends from the Earth surface down to a depth Z = 50 km. Its width is W = Z/sin α = 146 km. The coordinate x perpendicular to the line of strike, that is the distance from the trench over the Earth surface, is also shown. Positive and negative values identify points on the sides of the subducting (foot wall) and (hanging wall) plates, as indicated by the dash-dotted arrow. Figure 1. View largeDownload slide Schematicdepicting the 50 km thick outermost elastic layer, the 170 km thick asthenosphere and upper and lower mantle below. The locked portion of the fault, indicated by the dashed line, dips with angle α = 20° and extends from the Earth surface down to a depth Z = 50 km. Its width is W = Z/sin α = 146 km. The coordinate x perpendicular to the line of strike, that is the distance from the trench over the Earth surface, is also shown. Positive and negative values identify points on the sides of the subducting (foot wall) and (hanging wall) plates, as indicated by the dash-dotted arrow. For the sake of simplicity and with the aim of providing a first estimate of the geoid anomaly due to seismic and interseismic processes at subduction zones, we assume that thrust earthquakes occur at a constant recurrence period, T, and break uniformly the locked portion of the fault, with the same amount U of dip-slip at each earthquake. We thus write the stick-slip time history of the earthquake sequence as follows:   \begin{eqnarray} \delta u_e(t) = U\, \sum _{n=1}^{\infty } H(t-n\, T), \end{eqnarray} (1)where t is the time and H is the Heaviside function. Herein we assume that the first earthquake occurs at time t = T and, by convention, we start to count the seismic cycles from this time, This means that the nth seismic cycle is from n T to (n + 1) T. The initial period, from 0 to T, must be regarded as an arbitrary loading phase before the occurrence of the first earthquake. In addition to the stick-slip from each earthquake and according to Savage (1983), we consider also a back-slip distributed uniformly over the locked portion of the fault which grows linearly in time, starting from t = 0,   \begin{eqnarray} \delta u_b(t) = - \chi \, v\, t\, H(t) \end{eqnarray} (2)where χ is the seismic coupling coefficient and v is the convergence velocity between the subducting and overriding plates. This back-slip accounts for the strain accumulation process due to plate tectonics and constitutes the supplemental forcing that, once superimposed to the contribution from the steady state subduction (i.e. the long-term average motion that is not responsible for observable changes), describes the locking of the fault during the interseismic phase. Depending on the seismic coupling coefficient, this locking can be complete (χ = 1) or partial (χ < 1). For consistency, we equal the back-slip accumulated during one period, δub(t + T) − δub(t) = −χ v T to the stick-slip of one earthquake changed in sign, −U, and we obtain the following relation:   \begin{eqnarray} U = \chi \, v\, T. \end{eqnarray} (3) The dip-slip time history for which we will compute the geoid anomalies is the sum of the stick-slip and back-slip time histories:   \begin{eqnarray} \delta u(t) = \delta u_e(t)+ \delta u_b(t) = U\, f(t) \end{eqnarray} (4)with   \begin{eqnarray} f(t) = -\frac{t}{T}+\sum _{n=1}^{\infty } H(t-n\, T). \end{eqnarray} (5)In view of the arbitrariness of the choice of the initial time, the simulations that we show are significant only in the limit for t → ∞, after an infinite number of seismic cycles. We note that eq. (5) is equivalent to eq. (3) of Matsu’ura & Sato (1989) and, as in our case, it describes the slip history responsible for the periodic motion associated to the infinite sequence of seismic cycles within the seismogenic zone. Matsu’ura & Sato (1989) further account for the relative motion across the whole plate boundary between the subducting and overriding plates, their eq. (2). Nevertheless, as discussed by Matsu’ura & Sato (1989) and Sato & Matsu’ura (1993) for the case of old subduction zones, this relative motion is responsible only for a static contribution and, so, it does not yield any observable changes, in agreement with the assumptions underlying the back-slip approach proposed by Savage (1983). In this work, in order to focus on the fastest periodic changes within one seismic cycle, we do not consider the case of young subduction zones where this relative motion yields a growth rate of the topography of island arc-trench systems from one seismic cycle to the next one. 3 VISCOELASTIC RESPONSE From eq. (A6) and as discussed in Appendix A, within the framework of self-gravitating spherically symmetric viscoelastic Earth models with Maxwell rheology (Tanaka et al. 2006; Cambiotti et al. 2011c; Cambiotti & Sabadini 2015), the geoid anomaly at a point of the Earth surface, G, can be expressed as follows:   \begin{eqnarray} G(t) = (K\star \delta u)(t), \end{eqnarray} (6)where ⋆ stands for the time convolution and K is the viscoelastic Green function for a unitary, impulsive dip-slip over the whole locked portion of the fault. In order to compute the viscoelastic Green function, we consider a spherically symmetric Earth model based on PREM (table 3 of Dziewonski & Anderson (1981)) characterized by a 50 km thick elastic outermost layer and a 170 km thick (up to 220 km depth) asthenosphere, as portrayed in Fig. 1, with viscosity νa varied from 1017 to 1020 Pa s. The viscosities of the upper (from 220 km depth to the 670 km discontinuity) and lower (from the 670 km discontinuity to the core mantle boundary) mantle are set to 1021 and 3 × 1022 Pa s. This rheological stratification is consistent with most glacial isostatic adjustment and mantle convection studies (Cambiotti et al. 2011b; Tosi et al. 2005) and includes a low viscosity asthenosphere that mainly controls the relaxation of the earthquake-induced stress by viscous flow (Broerse et al. 2015). Furthermore, the initial density and bulk modulus are slightly modified in order to satisfy the Williamson–Adams equation (Wolf 1991; Wolf & Kaufmann 2000) and avoid Rayleigh–Taylor instabilities or compositional modes at large timescales (Plat & Jüttner 1995; Cambiotti & Sabadini 2010; Cambiotti et al. 2013). As shown by eqs (A12), (B10), (B11), (B13) and (A13a) and as discussed in Appendices A and B, eq. (6) with the dip-slip time history given by eq. (4) can be decomposed into three contributions:   \begin{eqnarray} G(t) = U\, \big (K_E\, f(t)+P(t)+C(t)\big ), \end{eqnarray} (7)where KE is the elastic Green function, P(t) is a periodic function with period T  \begin{eqnarray} P(t) = P(t+n\, T) \end{eqnarray} (8)such that P(n T) = 0, and C(t) has the following asymptotic limit   \begin{eqnarray} \lim _{t\rightarrow \infty }C(t) = C_\infty . \end{eqnarray} (9)The periodic function P(t) and this asymptotic limit can be expressed analytically, eqs (B14) and (B11), once the Laplace transform of the viscoelastic Green function has been evaluated at specific values of the Laplace variable. In view of the periodicity in time of the functions f(t) and P(t), the response of the Earth during the nth seismic cycle can be rewritten as   \begin{eqnarray} G_n(t) &= G(t+n\, T) = U \left(E(t)+P(t)+S_n(t)\right) & \end{eqnarray} (10)with t ∈ [0, T) and   \begin{eqnarray} E(t) = \left(1-t/T\right)\, K_E \end{eqnarray} (11a)  \begin{eqnarray} S_n(t) = C(t+n\, T)-K_E. & \end{eqnarray} (11b) We note that Sn(t) depends on the seismic cycle, while E(t) and P(t) do not. On the other hand, in view of eq. (9), we have   \begin{eqnarray} S_\infty = \lim _{n\rightarrow \infty } S_n(t) = C_\infty -K_E & \end{eqnarray} (12)and, then, the geoid anomaly after an infinite number of seismic cycles becomes a fully periodic function   \begin{eqnarray} G_\infty (t) = \lim _{n\rightarrow \infty }G_n(t)=U\, \big (E(t)+P(t)+S_\infty \big ). & \end{eqnarray} (13)From the latter formulation, we realize that E(t) and P(t) control the periodic changes of the geoid anomaly during one seismic cycle due to the instantaneous elastic response of the Earth and to the stress relaxation by viscous flow. Particularly, these periodic changes account for the instantaneous elastic response at the earthquake time, t = 0, and for the following changes during the interseismic phase which become zero just before the next earthquake, t = T−,   \begin{eqnarray} E(0)+P(0) = K_E \end{eqnarray} (14a)  \begin{eqnarray} E(T^-)+P(T^-) = 0. \end{eqnarray} (14b)We note also that the instantaneous elastic response, E(t), starts with the elastic response to stick-slip at the earthquake time and decreases linearly in time to zero just before the next earthquake. The constant S∞, instead, controls the static geoid anomaly over which these periodic changes are superimposed. This contribution results from stress relaxation by viscous flow due to the previous (infinite) seismic cycles and corresponds to the geoid anomaly evaluated just before the next earthquake that, as implied by eq. (12), reaches an asymptotic value after an infinite number of seismic cycles   \begin{eqnarray} G_\infty (T^{-}) = U\, S_\infty . \end{eqnarray} (15)This static geoid anomaly must not be confused with that caused by the steady state subduction. Our model, in fact, only focus on the strain accumulation and release process and does not account for the long-term relative motion between the subducting and overriding plates as in Matsu’ura & Sato (1989) and Sato & Matsu’ura (1993) and for the density heterogeneities of the subducting slab into the mantle and the induced mantle flow as in Ricard et al. (1993). Rather, as we are going to show in Section 5, it is a much smaller contribution resulting from the specific setting of the dip-slip time history describing the strain accumulation and release and that, although small and negligible, should be added to the static geoid anomaly caused by the steady state subduction. This analysis of the geoid anomaly after an infinite number of seismic cycles has allowed us to understand the physical meaning of E(t), P(t) and S∞. This understanding can now be used to gain insight into the response of the Earth after a finite number of seismic cycles. In this case, there is a net gain of the geoid anomaly accumulated during one interseismic period   \begin{eqnarray} G_{n+1}(t)-G_n(t) = S_{n+1}(t)-S_n(t) \end{eqnarray} (16)which involves only the function Sn(t). This net gain results from the fact that the stress relaxation by viscous flow due to all previous (finite) earthquakes and to the back-slip simulating the locking of the fault during the interseismic phase has not reached the steady state yet. Increasing the number of seismic cycles, in the limit for n → ∞, this net gain tends to zero and, then, there are no longer differences between a seismic cycle and the next one. 4 THE CASE OF A SIMPLE RELAXATION SPECTRUM In order to gain insight into the interaction between the stress relaxation and the dip-slip time history of the sequence of seismic cycles, it is useful to consider a heuristic Maxwell Earth model with a discrete relaxation spectrum. Within this assumption, the viscoelastic Green function can be expressed as follows:   \begin{eqnarray} K(t) = K_E\, \delta (t)+H(t)\, \sum _j k_j\, e^{-t/T_j} \end{eqnarray} (17)where δ(t) is the Dirac delta function, and kj and Tj are the residue and characteristic relaxation time of the jth relaxation mode. This specific form allows us to obtain the analytical expression for the time-dependent geoid anomaly. Particularly, as discussed in Appendix C, eq. (C5a), the periodic change due only to stress relaxation by viscous flow reads   \begin{eqnarray} P(t) = \sum _j r_j\, p_j(t) \end{eqnarray} (18)where rj = Tj Kj is the response to a stick-slip earthquake in the fluid limit, for t → ∞, associated to the jth relaxation mode (the so called strength) and pj(t) is a non-dimensional function describing the time evolution of this contribution during the seismic cycle   \begin{eqnarray} p_j(t) = \frac{1-e^{-t/T_j}}{1-e^{-T/T_j}}-\frac{t}{T}. \end{eqnarray} (19)This function is zero at the earthquake time, t = 0, increases to a maximum value $$p_j^*$$ at $$t=t^*_j$$, and then decreases to zero at the end of the seismic cycle, t = T. It is straightforward to obtain   \begin{eqnarray} p_j^* = \frac{1}{1-e^{-T/T_j}} +\frac{t_j^*-T_j}{T} \end{eqnarray} (20a)  \begin{eqnarray} t_j^* = T_j\, \left[\log \left(\frac{T}{T_j}\right)-\log \left(1-e^{-T/T_j}\right)\right]. \end{eqnarray} (20b) Fig. 2(a) shows the evolution in time of the function pj(t) for different relaxation times Tj = 0.01, 0.03, 0.1, 0.3, 1 and 3 T. We note that, increasing the relaxation times, the maximum amplitude $$p^*_j$$ decreases, while the time $$t_j^*$$ at which this maximum amplitude is attained increases. Fig. 2(b) shows also the maximum amplitude $$p_j^*$$ and the time $$t_j^*$$ at which this value is attained as function of the relaxation time Tj. We note that   \begin{eqnarray} \lim _{T_j/T\rightarrow 0} p_j^* = 1 \end{eqnarray} (21a)  \begin{eqnarray} \lim _{T_j/T\rightarrow \infty } p_j^* = 0. \end{eqnarray} (21b) This means that, for relaxation times smaller than the interseismic period, Tj ≪ T, the contribution of the relaxation mode corresponds almost to its strength, while for greater relaxation times, Tj ≫ T, its contribution is almost filtered out. This damping is already effective for relaxation times comparable with the interseismic period and, particularly, $$p_j^*\approx 1/2$$ for Tj ≈ T/5. As it concerns the time $$t_j^*$$ at which the contribution attains its maximum, it is always smaller than T/2 and this upper bound is reached only for relaxation times greater than the interseismic period, in the limit for Tj/T → ∞. In contrast, for relaxation time smaller than the interseismic period, Tj ≪ T, the time $$t_j^*$$ is about the relaxation time Tj scaled by log (T/Tj). Figure 2. View largeDownload slide (a) The contribution of the jth relaxation mode per unit strength during one seismic cycles of period, pj(t), for different relaxation times Tj = 0.01, 0.1, 0.3, 1 and 3 T. (b) The maximum amplitude of this contribution (left vertical axis, solid line) and the time at which this maximum is reached (right vertical axis, dashed line) as function of the relaxation time. Figure 2. View largeDownload slide (a) The contribution of the jth relaxation mode per unit strength during one seismic cycles of period, pj(t), for different relaxation times Tj = 0.01, 0.1, 0.3, 1 and 3 T. (b) The maximum amplitude of this contribution (left vertical axis, solid line) and the time at which this maximum is reached (right vertical axis, dashed line) as function of the relaxation time. 5 RESULTS We consider an interseismic period of T = 400 yr and, from eq. (3), assuming a partial seismic coupling of χ = 0.5 and a convergence velocity of v = 10 cm/yr−1, we set the stick-slip of each earthquake at U = 20 m. Considering a fault segment of (along strike) length L = 260 km and the shear modulus profile of PREM, each earthquake has moment magnitude Mw = 9.0. In the following we will focus first on the periodic geoid anomaly during one seismic cycle due to stress relaxation by viscous flow, U P(t), and we will consider the results for different asthenospheric viscosities in order to quantify the impact of the latter on the time-dependent pattern of this contribution. Then, once we will have clarified this issue, we will consider the whole periodic geoid anomaly, U (E(t) + P(t)), including the instantaneous elastic response, U E(t). Fig. 3 shows the periodic geoid anomaly during one seismic cycle due to stress relaxation by viscous flow, U P(t), along a transect over the Earth surface centred with the fault and perpendicular to the line of strike, according to Fig 1. This simulation is run for four asthenospheric viscosities (νa = 1017, 1018, 1019 and 1020 Pa s) and the geoid anomaly is shown as function of the distance from the trench, x, as depicted in Fig. 1. This periodic contribution is zero at the earthquake time, t = 0, and manifests itself only at later times according to the characteristic relaxation times of the Maxwell Earth model. For νa = 1017 Pa s, in the early phase of the seismic cycle, we have a positive geoid anomaly centred at around x = −150 km on the side of the overriding plate. This positive anomaly extends for about 600 ∼ 1000 km perpendicularly to the line of strike, mainly over the overriding plate. Elsewhere, instead, the geoid anomaly is negative, although small, less than a few millimetres at all times. The maximum geoid anomaly is about 1.7 cm already at t = 0.1 yr (about one month after the earthquake) and grows in time until t = 1 yr, when there is a maximum geoid anomaly of about 5.2 cm. At later times, the geoid anomaly begins to decrease and, at the same time, a local (positive) minimum close to the trench, at around x = −50 km on the side of the overriding plate, appears between two (positive) maxima, one on the side of the overriding plate, at around x = −200 km, and the other on the side of the subducting plate, at around x = 150 km. The main geoid anomaly remains positive in proximity of the trench until t = 10 yr and, successively, becomes negative in correspondence of the local minimum, reaching a minimum value of about −1.8 mm at y = 70 yr. Then, the amplitude of the geoid anomaly continues to decrease for all the distances from the trench until it is zero everywhere at the end of the seismic cycle. Figure 3. View largeDownload slide Periodic geoid anomaly during one seismic cycle due to stress relaxation by viscous flow, U P(t), along a transect over the Earth surface centred with the fault and perpendicular to line of strike. The interseismic period and the stick-slip of each earthquake are set to T = 400 yr and U = 20 m. The locked fault crosses the whole elastic layer, down to a depth Z = 50 km. Each panel shows the geoid anomaly as function of the distance from the trench, x, at increasing time (from top to bottom and from left to right; see the label in each panel) for asthenosphere viscosities of 1017, 1018, 1019 and 1020 Pa s (solid, dashed, dashed-dotted and dotted lines, respectively). The vertical grey line indicates the trench, at x = 0. Figure 3. View largeDownload slide Periodic geoid anomaly during one seismic cycle due to stress relaxation by viscous flow, U P(t), along a transect over the Earth surface centred with the fault and perpendicular to line of strike. The interseismic period and the stick-slip of each earthquake are set to T = 400 yr and U = 20 m. The locked fault crosses the whole elastic layer, down to a depth Z = 50 km. Each panel shows the geoid anomaly as function of the distance from the trench, x, at increasing time (from top to bottom and from left to right; see the label in each panel) for asthenosphere viscosities of 1017, 1018, 1019 and 1020 Pa s (solid, dashed, dashed-dotted and dotted lines, respectively). The vertical grey line indicates the trench, at x = 0. The pattern of the geoid anomaly for νa = 1018 Pa s is similar, although delayed in time by about a factor of 10 with respect to the case of νa = 1017 Pa s. This means that the relaxation times of the Maxwell Earth model are mainly controlled by the asthenospheric viscosity. The maximum geoid anomaly is about 1.7 cm after one year, t = 1 yr, and grows up to 5.2 cm at 10 yr. Afterwards, the local positive minimum close to the trench and between the two local maxima begins to appear. Then, it becomes negative at 150 yr and reaches its minimum value of −0.3 cm at 250 yr, when the amplitude of the geoid anomaly is already decreasing for all distance from the trench until it is zero everywhere at the end of the seismic cycle. The cases of νa = 1019 and 1020 Pa s are simpler because the geoid anomalies consist mainly of positive geoid anomalies in proximity of the trench, the amplitudes of which grow in time until t = 70 and 150 yr, when they reach the maximum values of 4.2 and 0.9 cm, respectively. At later times, the amplitudes of these positive geoid anomalies begin to decrease for all distances from the trench and are zero at the end of the seismic cycle. It is noteworthy that the maximum amplitudes of the geoid anomaly over the whole seismic cycle for both cases of νa = 1017 and 1018 Pa s are similar (about 5.2 cm), although they are attained at different times (t = 1 and 10 yr, respectively). In contrast, the maximum amplitude of the geoid anomaly obtained for νa = 1019 Pa s is less by about the 20 and 80 per cent than this value, only 4.2 and 0.9 cm at t = 70 and 150 yr, respectively. This reduction of the amplitude of the geoid anomaly by increasing the asthenospheric viscosity is consistent with the discussion of Section 4 and is explained considering that the main relaxation modes of the Maxwell Earth models with νa = 1017 and 1018 Pa s have relaxation modes with characteristic relaxation times smaller than the interseismic period T = 400 yr and, then, they contribute with almost their complete strengths during the interseismic phase. In contrast, the main relaxation modes of the Maxwell Earth model with νa = 1019 have larger characteristic relaxation times and their contribution is partially filtered out. This is further confirmed by the simulation obtained for asthenospheric viscosity of 1020 Pa s. In this case, the temporal evolution is similar to that obtained for νa = 1019 Pa s, but the maximum amplitude is even smaller. Having discussed the periodic geoid anomaly due to stress relaxation by viscous flow, we now consider the whole periodic geoid anomaly, U (E(t) + P(t)), including the instantaneous elastic response as shown in Fig. 4. The geoid anomaly at the earthquake time, t = 0, is the same for all asthenospheric viscosities and corresponds to the elastic response to the stick-slip earthquake, U KE. It is characterized by a positive geoid anomaly in proximity of the trench, with a maximum of about 3.8 cm on the side of the overriding plate, at around x = −40 km, and a negative geoid anomaly on the side of the overriding plate, with a minimum of about −1.3 cm at around x = −180 km. According to eq. (11a), this elastic response decreases linearly in time during the interseismic phase and yields zero at the end of the seismic cycle. Figure 4. View largeDownload slide As Fig. 3, but including the instantaneous elastic response, U (E(t) + P(t)). Figure 4. View largeDownload slide As Fig. 3, but including the instantaneous elastic response, U (E(t) + P(t)). At the early stage, the stress relaxation increases the positive geoid anomaly caused by the elastic response in proximity of the trench, and, being wider in space, compensates the negative geoid anomaly at around x = −150 km, on the side of the overriding plate. In particular, this compensation of the negative geoid anomaly due to the instantaneous elastic response occurs at t = 0.1, 1, 10 and 150 yr for νa = 1017, 1018, 1019 and 1020 Pa s. The maximum geoid anomaly attained during the seismic cycle is 6.6 cm at t = 1 yr, 6.5 cm at 10 yr and 5.6 cm at 40 yr for νa = 1017, 1018 and 1019 Pa s. For the case of νa = 1020, instead, the maximum geoid anomaly decreases at later times due to the fact that the positive contribution due to stress relaxation by viscous flow, U P(t), is smaller than the reduction, linearly in time, of the elastic response to the stick- and back-slips, U E(t). The short wavelengths of the elastic response remain visible at all times and make the maximum geoid anomaly always greater than that from the periodic contribution due to stress relaxation by viscous flow. This pattern is not altered even when the local minimum close to the trench between the two maxima appears at t = 2 and 20 yr for νa = 1017 and 1018 Pa s (Fig. 3). Furthermore, we note that the periodic geoid anomaly for νa = 1017 Pa s becomes smaller than that for νa = 1018, 1019 and 1020 Pa s between 2 to 4 yr, 4 to 10 yr and 10 to 20 yr. Particularly, at t = 200 yr which is one half the interseismic period, T/2, the amplitude of the geoid anomaly for νa = 1017 Pa s does not exceed 0.7 cm along the whole transect. At this time, instead, the geoid anomalies for νa = 1018, 1019 and 1020 Pa s reach maximum values of 1.6, 3.2 and 2.4 cm close to the trench. This results from the faster decay of the geoid anomaly due to stress relaxation for the lowest asthenospheric viscosity of 1017 Pa s. At later times, the geoid anomalies for all the asthenospheric viscosities decay to zero everywhere, as expected on the basis of the results obtained in Section 4 for the heuristic case of a discrete relaxation spectrum. In fact, the contribution from each relaxation mode reaches its maximum amplitude before one half the interseismic period, T/2 = 200 yr. In order to provide a clearer description of our results, we show in Fig. 5 the periodic geoid anomaly as function of time at two representative points on the side of the overriding plates, at x = −150 and −50 km. The point at x = −150 km is about where the maximum geoid anomaly during the early stage of the seismic cycle is located. The point at x = −50 km, instead, is about where the local minimum between the two maxima for the cases of νa = 1017 and 1018 Pa s is located. As expected, the periodic geoid anomaly due to only stress relaxation resembles that of the functions pj controlling the time evolution of the contributions of the relaxation modes, within the framework of the heuristic Earth model with a discrete relaxation spectrum. Differences, like the change in sign in Fig. 5(b) for νa = 1017 Pa s and in Fig. 5(d) for νa = 1017 and 1018 Pa s and the different concavities after the early phase of the seismic cycle for νa = 1017 and 1018 Pa s, are the result of the superimposition of relaxation modes with strengths of opposite signs. Figure 5. View largeDownload slide (a,c) Periodic geoid anomaly during one seismic cycle due to both the instantaneous elastic response and stress relaxation by viscous flow U (E(t) + P(t)); (b,d) as in (a,c) but limited to relaxation by viscous flow U P(t), at the distances from the trench x = −150 (a,b) and −50 km for (c,d). Solid, dashed, dashed-dotted and dotted lines refer to the case of asthenospheric viscosities of νa = 1017, 1018, 1019 and 1020 Pa s. The thick solid line in panels (a) and (c) indicates the instantaneous elastic response, U E(t). Figure 5. View largeDownload slide (a,c) Periodic geoid anomaly during one seismic cycle due to both the instantaneous elastic response and stress relaxation by viscous flow U (E(t) + P(t)); (b,d) as in (a,c) but limited to relaxation by viscous flow U P(t), at the distances from the trench x = −150 (a,b) and −50 km for (c,d). Solid, dashed, dashed-dotted and dotted lines refer to the case of asthenospheric viscosities of νa = 1017, 1018, 1019 and 1020 Pa s. The thick solid line in panels (a) and (c) indicates the instantaneous elastic response, U E(t). At the end, for the sake of completeness, we show in Fig. 6 also the static geoid anomaly from the steady-state response after the infinite series of seismic cycles, U S∞, on which the periodic change shown in Fig. 4 is superimposed. As anticipated at the end of Section 3, this contribution is smaller, with amplitudes of at most 3.2 cm, than the static geoid anomaly due to the long term average steady state motion at subduction zones and to the mantle density heterogeneities and the induced mantle flow, which is of the order of tens of meters (Matsu’ura & Sato 1989; Sato & Matsu’ura 1993; Zhong & Gurnis 1992; Marotta et al. 2006). Figure 6. View largeDownload slide The static contribution due to stress relaxation by viscous flow after an infinite number of seismic cycles, U S∞, for asthenospheric viscosities of νa = 1017, 1018, 1019 and 1020 Pa s (solid, dashed, dashed-dotted, and dotted lines). Figure 6. View largeDownload slide The static contribution due to stress relaxation by viscous flow after an infinite number of seismic cycles, U S∞, for asthenospheric viscosities of νa = 1017, 1018, 1019 and 1020 Pa s (solid, dashed, dashed-dotted, and dotted lines). 6 CONCLUSION We have extended the approach proposed by Savage (1983) from gravitating and flat to self-gravitating and spherically symmetric Earth models in order to quantify the geoid anomaly caused by the strain accumulation and release at subduction zones. Different from Savage (1983) that considered periodic topography changes during the seismic cycles, we have focused on the geoid anomaly and obtained analytical expressions for the periodic geoid anomaly after an infinite number of seismic cycles, as well as for the a smaller static contribution from the steady-state response over which this periodic change is superimposed. These analytical expressions do not require to calculate the viscoelastic response to each earthquake in the time domain, as done instead by Savage (1983) on the basis of the solution provided by Thatcher & Rundle (1981). Instead, they are based on the Correspondence Principle between the viscoelastic problem in the time domain and an equivalent elastic problem in the Laplace domain, and the analysis of the singularities of the geoid anomaly in the Laplace domain discussed in Appendix B. Further, within the heuristic assumption of a Maxwell Earth model with a discrete relaxation spectrum, we have understood how the contribution to the periodic geoid anomaly from relaxation modes with characteristic relaxation times comparable or greater than the interseismic period is damped, until it is completely filtered out when the relaxation times are greater than the interseismic period. We have applied this methodology to the case of an infinite sequence of megathrust earthquakes with moment magnitude Mw = 9 and an interseismic period of T = 400 yr, and computed the periodic geoid anomaly during the seismic cycle. The elastic response at the earthquake time yields a positive geoid anomaly with a maximum value of 3.8 cm close to the trench and on the side of the overriding plate and a negative geoid anomaly further away from the trench with a minimum value of −1.3 cm (Fig. 4 at t = 0 yr). At later times, this instantaneous elastic response decreases to zero linearly in time due to the elastic contribution associated with the locking of the fault during the interseismic phase. During the interseismic phase, the periodic geoid anomaly is affected also by stress relaxation by viscous flow. The latter provides a peculiar pattern of global upwarping of the geoid that is broader in space compared to the elastic contribution and has its maximum close to the trench and on the side of the overriding plate (Fig. 3). This positive geoid anomaly grows in the early stage up to values that can be two times greater than the maximum amplitude of the elastic response to the earthquakes for low asthenospheric viscosities. Particularly, the maximum geoid anomalies due to stress relaxation attain maximum values of 5.2 at 1 and 10 yr after the earthquake for asthenospheric viscosities of 1017 and 1018 Pa s. In contrast, for the case of an asthenospheric viscosities of 1019 and 1020 Pa s, the maximum positive geoid anomalies attain maximum values of 4.2 and 0.9 cm at 70 and 150 yr after the earthquake. This is consistent with the finding that the interaction between stress relaxation by viscous flow and the periodic dip-slip history describing strain accumulation and release damps the contribution from those relaxation modes which have relaxation times comparable or greater than the interseismic period (Fig. 2). Once the geoid anomaly due to stress relaxation has attained its maximum, it starts to decrease to zero. In this second phase of the seismic cycle, we note that the lowermost asthenospheric viscosity of 1017 Pa s yields the fastest collapse of the geoid maximum which is splitted into two reduced maxima intermingled by a local minimum close to the trench, finally developing into a negative geoid anomaly starting from t = 10 yr after the earthquake. These results constitutes a first step towards the interpretation of time-dependent geoid anomalies in the surroundings of subduction zones and can be used as a theoretical framework for explaining data time series from the GRACE and GOCE satellite missions even at those locations where no megathurst earthquakes have occurred during the time in which these two missions have been operative, from 2002 to present and from 2009 to 2013. Particularly, an improved time series analysis by Meyer et al. (2016), based on the optimization of the GRACE data, leads to the new release AIUB (Astronomical Institute of the University of Bern) — RL02 and allows to estimate the standard deviation (formal errors) for the time-dependent geoid evaluation up to the harmonic degree ℓ = 60. This standard deviation, evaluated from the difference between the time-dependent geoid and a static reference geoid, which is concordant with our approach where the time-dependent geoid must be considered as superimposed to the static geoid due to steady-state subduction, is about 1 mm (fig. 2 of Meyer et al. 2016) and indicates that our co-, post- and interseismic signals are detectable from GRACE, since in the interseismic phase the geoid anomaly grows to several centimetres on timescales ranging from one year to a few decades, depending on the asthenospheric viscosities. Although our results take into account all the wavelengths of the geoid anomaly, in view of its large spatial scale, the filtering in the bandwidth up to the harmonic degree ℓ = 60 of GRACE does not affect this conclusion. Indeed, this filtering reduces the amplitude of the geoid anomalies due to the elastic response, U E, and to stress relaxation by viscous flow, U P, by about a factor of 4 and 2, respectively, with the largest reduction affecting the elastic response because of its shorter wavelengths. The comparison of our modelling with space gravity data, in addition to proper wavelength filtering, requires specific assumptions on the locked portion of the fault during the interseismic phase, which does not necessarily coincide with that adopted in this work (extending form the Earth surface down to the outermost 50 km), and the knowledge of past seismicity back in time as far as the post-seismic contribution significantly affects the data time series. The sensitivity analysis to different locking of the fault and to the realistic earthquake sequences will be the subject of future works. Acknowledgements This research was partially supported by the CAS/CAFEA international partnership program for creative research teams (No. KZZDEW-TZ-19). DY acknowledges support from Geochemistry and CISE from National Science Foundation. REFERENCES Broerse T., Riva R., Vermeersen B., 2011. Ocean contribution to seismic gravity changes: the sea level equation for seismic perturbations revisited, Geophys. J. Int. , 199, 1094– 1109. Google Scholar CrossRef Search ADS   Broerse T., Riva R., Simons W., Govers R., Vermeersen B., 2015. Postseismic grace and GPS observations indicate a rheology contrast above and below the Sumatra slab, J. geophys. Res. , 120, 5343– 5361. Google Scholar CrossRef Search ADS   Cambiotti G., Sabadini R., 2010. The compressional and compositional stratifications in Maxwell Earth models: the gravitational overturning and the long-period tangential flux, Geophys. J. Int. , 180, 475– 500. Google Scholar CrossRef Search ADS   Cambiotti G., Sabadini R., 2013. Gravitational seismology retrieving centroid moment tensor solution of the 2011 Tohoku earthquake, J. geophys. Res. , 118, 183– 194. Google Scholar CrossRef Search ADS   Cambiotti G., Sabadini R., 2015. On the seismic perturbation due to a fault system: its evaluation beyond the epicentral reference frame, Geophys. J. Int. , 203, 943– 959. Google Scholar CrossRef Search ADS   Cambiotti G., Bordoni A., Sabadini R., Colli L., 2011a. Grace gravity data help constraining seismic models of the 2004 Sumatran earthquake, J. geophys. Res. , 116, B10403. Google Scholar CrossRef Search ADS   Cambiotti G., Ricard Y., Sabadini R., 2011b. New insights into mantle convection true polar wander and rotational bulge readjustment, Earth planet. Sci. Lett. , 310, 538– 543. Google Scholar CrossRef Search ADS   Cambiotti G. Y. R., Sabadini R., 2011c. Ice age True PolarWander in a compressible and non-hydrostatic Earth, Geophys. J. Int. , 183, 1248– 1264. Google Scholar CrossRef Search ADS   Cambiotti G., Klemann V., Sabadini R., 2013. Compressible viscoelastodynamics of a spherical body at long timescales and its isostatic equilibrium, Geophys. J. Int. , 193, 1071– 1082. Google Scholar CrossRef Search ADS   Chen J., Wilson C., Tapley B., Grand S., 2007. Grace detects coseismic and postseismic deformation from the sumatra-andaman earthquake, Geophys. Res. Lett. , 34, L13302. De Linage C., Rivera L., Hinderer J., Boy J.-P., Rogister Y., Lambotte S., Biancale R., 2009. Separation of coseismic and postseismic gravity changes for the 2004 Sumatra–Andaman earthquake from 4.6 yr of GRACE observations and modelling of the coseismic change by normal-modes summation, Geophys. J. Int. , 176, 695– 714. Google Scholar CrossRef Search ADS   Dziewonski A., Anderson D., 1981. Preliminary Reference Earth Model, Phys. Earth planet. Inter. , 25, 297– 356. Google Scholar CrossRef Search ADS   Fuchs M., Hooper A., Broerse T., Bouman J., 2016. Distributed fault slip model for the 2011 Tohoku-Oki earthquake from GNSS and GRACE/GOCE satellite gravimetry, J. geophys. Res. , 121, 1114– 1130. Google Scholar CrossRef Search ADS   Gross R., Chao B., 2001. The gravitational signature of earthquakes, in Gravity, Geoid and Geodynamics 2000 , 123, pp. 205– 210, ed. Sideris M.G., International Association of Geodesy Symposia. Google Scholar CrossRef Search ADS   Han S., Sauber J., Luthcke S., Ji C., Pollitz F., 2008. Implications of postseismic gravity change following the great 2004 Sumatra–Andaman earthquake from the regional harmonic analysis of grace intersatellite tracking data, J. geophys. Res. , 113, B11413. Google Scholar CrossRef Search ADS   Marotta A., Spelta E., Rizzetto C., 2006. Gravity signature of crustal subduction inferred from numerical modelling, Geophys. J. Int. , 166, 923– 938. Google Scholar CrossRef Search ADS   Matsuo K., Heki K., 2011. Coseismic gravity changes of the 2011 Tohoku-Oki earthquake from satellite gravimetry, Geophys. Res. Lett. , 38, L00G12. Google Scholar CrossRef Search ADS   Matsu’ura M., Sato T., 1989. A dislocation model for the earthquake cycle at convergent plate boundaries, Geophys. J. Int. , 96, 23– 32. Google Scholar CrossRef Search ADS   Meyer U., Jäggi A., Jean Y., Beutler G., 2016. Aiub-rl02: an improved time-series of monthly gravity fields from grace data, Geophys. J. Int. , 205, 1196– 1207. Google Scholar CrossRef Search ADS   Plat H.-P., Jüttner H.-U., 1995. Rayleigh-Taylor instabilities of a self-gravitating Earth, J. Geodyn. , 20, 267– 288. Google Scholar CrossRef Search ADS   Pollitz F., Bürgmann R., Banerjee P., 2006. Post-seismic relaxation following the great 2004 Sumatra–Andaman earthquake on a compressible self-gravitating earth, Geophys. J. Int. , 167, 397– 420. Google Scholar CrossRef Search ADS   Ricard Y., Richards M., Lithgow-Bertelloni C., Lestunff Y., 1993. A geodynamic model of mantle density heterogeneity, J. geophys. Res. , 98, 21 895– 21 909. Google Scholar CrossRef Search ADS   Sato T., Matsu’ura M., 1993. A kinematic model for evolution of island arc-trench systems, Geophys. J. Int. , 114, 512– 530. Google Scholar CrossRef Search ADS   Savage J.C., 1983. A dislocation mode of strain accumulation and release at a subduction zone, J. geophys. Res. , 88, 4984– 4996. Google Scholar CrossRef Search ADS   Savage J.C., Burford R.O., 1973. Geodetic determination of relative plate motion in Central California, J. geophys. Res. , 78, 832– 845. Google Scholar CrossRef Search ADS   Segall P., 2010. Earthquake and Volcano Deformation , Princeton University Press. Google Scholar CrossRef Search ADS   Song T., Simons M., 2003. Large Trench-Parallel Gravity Variations Predict Seismogenic Behavior in Subduction Zones, Science , 301, 630– 633. Google Scholar CrossRef Search ADS PubMed  Sun W., Okubo S., 2004. Co-seismic deformations detectable by satellite gravity missions, J. geophys. Res. , 109, B04405. Tanaka Y., Okuno J., Okubo S., 2006. A new method for the computation of global viscoelastic post-seismic deformation in a realistic earth model (I)–vertical displacement and gravity variation, Geophys. J. Int. , 164, 273– 289. Google Scholar CrossRef Search ADS   Tanaka Y., Hasegawa T., Tsuruoka H., Klemann V., Martinec Z., 2015. Spectral-finite element approach to post-seismic relaxation in a spherical compressible Earth: application to gravity changes due to the 2004 Sumatra–Andaman earthquake, Geophys. J. Int. , 200, 299– 321. Google Scholar CrossRef Search ADS   Thatcher W., Rundle J., 1979. A model for the earthquake cycle in underthrust zones, J. geophys. Res. , 84, 5540– 5556. Google Scholar CrossRef Search ADS   Thatcher W., Rundle J., 1981. A viscoelastic model for periodically recurring earthquakes in subduction zones (abstract), EOS, Trans. Am. geophys. Un. , 62, 1031. Tosi N., Sabadini R., Marotta A. M., Vermeersen L., 2005. Simultaneous inversion for the Earth’s mantle viscosity and ice mass imbalance in Antarctica and Greenland, J. geophys. Res. , 110, B07402. Google Scholar CrossRef Search ADS   Wang L., Shum C.F.J.S., Tapley B., Dai C., 2012. Coseismic and postseismic deformation of the 2011 Tohoku-Oki earthquake constrained by GRACE gravimetry, Geophys. Res. Lett. , 39, L07301. Wolf D., 1991. Visco-elastodynamics of a stratified compressible planet: incremental field equations and short- and long-time asymptotes, Geophys. J. Int. , 104, 401– 417. Google Scholar CrossRef Search ADS   Wolf D., Kaufmann G., 2000. Effects due to compressional and compositional density stratification on load-induced Maxwell-viscoelastic perturbations, Geophys. J. Int. , 140, 51– 62. Google Scholar CrossRef Search ADS   Zhong S., Gurnis M., 1992. Viscous flow model of a subduction zone with a faulted lithosphere: Long and short wavelength topography, gravity, and geoid, Geophys. Res. Lett. , 19, 1891– 1894. Google Scholar CrossRef Search ADS   APPENDIX A: ELASTIC AND VISCOELASTIC GREEN’s FUNCTIONS Within the framework of self-gravitating spherically symmetric elastic Earth models and adapting the formalism presented in Cambiotti & Sabadini (2015) for the case of a time-dependent uniform slip over the fault surface (i.e. the slip varies through time but it does not depend on the position over the fault), the geoid anomaly at a point of the Earth surface can be expressed as follows:   \begin{eqnarray} G(\theta ,\varphi ,t) = \sum _{j=1}^{2} \delta u_j(t)\, K_E^j(\theta ,\varphi ) \end{eqnarray} (A1)where δuj are the along strike (j = 1) and along dip (j = 2) components of the uniform slip, and $$K^j_E(\theta ,\varphi )$$ are the elastic Green functions for the geoid anomaly at colatitude θ and longitude φ given by   \begin{eqnarray} K^j_E(\theta ,\varphi )=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell } K_{E\ell m}^j\, Y_{\ell m}(\theta ,\varphi ). \end{eqnarray} (A2)Here, Yℓm are the normalized real spherical harmonics of degree ℓ and order m, and $$K_{E\ell m}^j$$ are the spherical harmonic coefficients of $$K^j_E(\theta ,\varphi )$$. The latter can be obtained as follows:   \begin{eqnarray} K_{E\ell m}^j = a^{-2}\, \sum _{p=0}^{2}\int _{\mathcal {S}} k_{\ell ,p}(a,r^{\prime })\, \Delta _{\ell m,pj}^D(\mathbf {r}^{\prime })\, \mathrm{d}S^{\prime } \end{eqnarray} (A3)where a is the Earth radius, $$\mathcal {S}$$ is the fault surface, dS΄ is the infinitesimal area, r΄ is the vector position identifying a point of the fault surface, kℓ, p(a, r΄) is the (elastic) gravitational seismic Love number governing the polar (p = 0), bipolar (p = 1) and quadrupolar (p = 2) patterns of the response at the Earth surface to a point-like seismic force at the radial distance from the Earth centre r΄ = ‖r΄‖, and Δℓm, pj is the following factor:   \begin{eqnarray} \Delta _{\ell m,pj}(\mathbf {r}^{\prime })= \frac{D_j^p(\mathbf {r}^{\prime })\, F_{\ell m}^p(\theta ^{\prime },\varphi ^{\prime })+D_j^{-p}(\mathbf {r}^{\prime })\, F_{\ell m}^{-p}(\theta ^{\prime },\varphi ^{\prime })}{1+\delta _{p0}} \end{eqnarray} (A4)which accounts for the geometry (dip and strike angles) and location (colatitude θ΄ and longitude φ΄) of the infinitesimal surface element, with $$F_{\ell m}^p$$ and $$D_{j}^p$$ given by eqs (31) and (33) of Cambiotti & Sabadini (2015). In order to perform the integration over depth z΄ (or the radial distance from the Earth surface r΄ = a − z΄) we interpolate the seismic Love numbers sampled every 1 km, but for the shallowest 3 km where the sampling is gradually refined up to 0.1 km. In the main text we will not explicit the dependence on the colatitude and longitude of the point at which the geoid anomaly is computed and, as far as we will focus only on dip-slip faults, we will omit the indexing for distinguishing between strike-slip and dip-slip and we simply write   \begin{eqnarray} G(t) = K_E\, \delta u(t) \end{eqnarray} (A5)where it is understood that δu and KE correspond to δu2 and $$K^2_E$$ in the formulation used in this appendix. According to the correspondence principle (Tanaka et al. 2006), these results can be generalized to Maxwell Earth models as follows:   \begin{eqnarray} G(t) = (K\star \delta u)(t) \end{eqnarray} (A6)denoting the time convolution between the slip time history δu and the viscoelastic Green function K. The latter is defined as the inverse Laplace transform of the equivalent elastic Green functions $$\tilde{K}$$ in the Laplace domain   \begin{eqnarray} K(t)=\mathcal {L}^{-1}[\tilde{K}(s)] \end{eqnarray} (A7)where $$\tilde{K}(s)$$ is obtained solving the equivalent elastic problem where the shear modulus μ is replaced by the equivalent shear modulus $$\hat{\mu }(s),$$  \begin{eqnarray} \hat{\mu }(s) = \frac{\mu s}{s+1/\tau }, \end{eqnarray} (A8)providing the dependence on the Laplace variable s. The lambda parameter λ, instead, is replaced by   \begin{eqnarray} \hat{\lambda }(s) = \kappa - \frac{2}{3}\, \hat{\mu }(s). \end{eqnarray} (A9) In order to distinguish between the instantaneous elastic response and the contributions caused by stress relaxation by viscous flow, the viscoelastic Green function is usually decomposed as follows:   \begin{eqnarray} K(t) = K_E\, \delta (t)+K_V(t) \end{eqnarray} (A10)with   \begin{eqnarray} K_V(t)=\mathcal {L}^{-1}[\tilde{K}_V(s)] \end{eqnarray} (A11a)  \begin{eqnarray} \tilde{K}_V(s) = \tilde{K}(s)-K_E. \end{eqnarray} (A11b) Then, substituting eq. (A10) into eq. (A6), we obtain   \begin{eqnarray} G(t) = G_E(t)+G_V(t), \end{eqnarray} (A12)where GE and Gv are the instantaneous elastic response and the response due to only stress relaxation by viscous flow,   \begin{eqnarray} G_E(t) =\delta u(t)\, K_E \end{eqnarray} (A13a)  \begin{eqnarray} G_V(t) =\mathcal {L}^{-1}[\tilde{G}_V(s)] \end{eqnarray} (A13b)with   \begin{eqnarray} \tilde{G}_V(s)=\tilde{K}_V(s)\, \delta \tilde{u}(s). \end{eqnarray} (A14)Here, we have made use of the fact that the time convolution between two functions corresponds to the product of their Laplace transforms in the Laplace domain. In the summation over the harmonic degrees in eq. (A2) we have considered the harmonic degrees up to 5 × 105 and 2 × 103 in order to obtain the elastic, GE, and viscous, Gv, Green functions, respectively. We note that the viscous Green function requires a much lower truncation than the elastic one because it reflects the stress relaxation by viscous flow that occurs, within the assumption of the present work, at depths greater than 50 km. APPENDIX B: GEOID ANOMALY In order to estimate the geoid anomaly due to only stress relaxation by viscous flow, we first consider the Laplace transform of the dip-slip time history that, by making use of the expression for geometric series, we rearrange as follows:   \begin{eqnarray} \delta \tilde{u}(s) = U\, \left(-\frac{1}{s^2\, T}+\sum _{n=1}^{\infty } \frac{e^{-n\, s\, T}}{s}\right)=U\, \frac{\tilde{g}(s)}{s} \end{eqnarray} (B1)with   \begin{eqnarray} \tilde{g}(s) = -\frac{1}{1-e^{s\, T}}-\frac{1}{s\, T} \end{eqnarray} (B2)Then, substituting eq. (B1) into eq. (A14), we obtain   \begin{eqnarray} \tilde{G}_V(s) = U\, \frac{\tilde{K}_V(s)\, \tilde{g}(s)}{s}. \end{eqnarray} (B3) In order to perform the inverse Laplace transform of $$\tilde{G}_V(s)$$, we must first investigate its singularities in the Laplace domain. They are depicted in Fig. B1. The function $$\tilde{K}_V(s)$$ has singularities on the negative half of the real axis of the Laplace domain and, particularly, they are confined to the interval $$(-1/\tau _{\rm {min}})\subset \mathbb {R}^-$$. These singularities depend on the specific viscoelastic stratification of the Earth model and can be isolated poles, clusters of poles and continuous not-Lipschitzian intervals (Tanaka et al. 2006; Cambiotti et al. 2011c). The function $$\tilde{g}(s)$$, instead, has singularities at the isolated poles s = i ωn (n = ±1, …, ±∞) with   \begin{eqnarray} \omega _n = \frac{2\, \pi \, n}{T} \end{eqnarray} (B4)and its residues are   \begin{eqnarray} \lim _{s\rightarrow \infty } \tilde{g}(s)\, (s-i\, \omega _n) = \frac{1}{T}. \end{eqnarray} (B5)We note that s = 0 is not a singularity for $$\tilde{g}(s)$$ and, in fact, its limit for s → 0 exists   \begin{eqnarray} \lim _{s\rightarrow 0} \tilde{g}(s) = -\frac{1}{2}. \end{eqnarray} (B6) Figure B1. View largeDownload slide Schematic depicting the singularities of the function $$\tilde{K}_V(s)\, \tilde{g}(s)$$ in the Laplace domain. The black circles indicate the poles along the imaginary axis due to $$\tilde{g}(s)$$, which are at s = i ωn (n = ±1, …, ±∞). The white circles and rectangles indicate the poles and the not-Lipschitzian intervals due to $$\tilde{K}_V(s)$$. The latter lay on the negative half of the real axis of the Laplace domain and, particularly, within the interval $$(-1/\tau _{\rm {min}},0)\subset \mathbb {R}^-$$ identified by the thick solid line. The dashed rectangle is the closed contour Γ that includes the origin of the Laplace domain and all the singularities due to $$\tilde{K}_V(s)$$. Figure B1. View largeDownload slide Schematic depicting the singularities of the function $$\tilde{K}_V(s)\, \tilde{g}(s)$$ in the Laplace domain. The black circles indicate the poles along the imaginary axis due to $$\tilde{g}(s)$$, which are at s = i ωn (n = ±1, …, ±∞). The white circles and rectangles indicate the poles and the not-Lipschitzian intervals due to $$\tilde{K}_V(s)$$. The latter lay on the negative half of the real axis of the Laplace domain and, particularly, within the interval $$(-1/\tau _{\rm {min}},0)\subset \mathbb {R}^-$$ identified by the thick solid line. The dashed rectangle is the closed contour Γ that includes the origin of the Laplace domain and all the singularities due to $$\tilde{K}_V(s)$$. We can take advantage of the fact that we have analytical expressions for the poles and the residues of $$\tilde{g}(s)$$ and that $$\tilde{K}_V(s)$$ is analytical at these poles decomposing the product $$\tilde{K}_V(s)\, \tilde{g}(s)$$ as follows:   \begin{eqnarray} \tilde{K}_V(s)\, \tilde{g}(s) = \tilde{p}(s)+\tilde{c}(s), \end{eqnarray} (B7)where $$\tilde{p}(s)$$ accounts for all the poles at s = i ωn:   \begin{eqnarray} \tilde{p}(s) = \frac{1}{T}\, \sum _{n=1}^\infty \left(\frac{\tilde{K}_V(i\, \omega _n)}{s-i\, \omega _n}+ \frac{\tilde{K}_V^*(i\, \omega _{n})}{s+i\, \omega _n}\right). \end{eqnarray} (B8)Here * stands for the complex conjugate and we have made use of the fact that the Laplace transforms of real functions has the following property: $$K_V(s^*) = K_V^*(s)$$. The function $$\tilde{c}(s)$$, instead, is simply the difference between the product $$\tilde{K}^V(s)\, \tilde{g}(s)$$ and $$\tilde{p}(s),$$  \begin{eqnarray} \tilde{c}(s) = \tilde{K}_V(s)\, \tilde{g}(s)-\tilde{p}(s). \end{eqnarray} (B9)This function has the same singularities of $$\tilde{K}_V(s)$$, but it has not those from $$\tilde{g}(s)$$ because they are removed by subtracting $$\tilde{P}(s)$$. Indeed, the residues of $$\tilde{C}(s)$$ at s = i ωn are zero. Within this decomposition, the geoid anomaly due to stress relaxation by viscous flow, eq. (A13b), that is the inverse Laplace transform of eq. (B3), can be written as follows:   \begin{eqnarray} G_V(t) = U \left(P(t)+C(t)\right)\!, \end{eqnarray} (B10)where P(t) is the inverse Laplace transform of $$\tilde{p}(s)/s$$, which reads   \begin{eqnarray} P(t) \!=\! \mathcal {L}^{-1}\big [\tilde{p}(s)/s\big ] \!=\! \sum _{n=1}^\infty \big [I_n\, (1\!-\!\cos (\omega _n\, t))\!+\!R_n\, \sin (\omega _n\, t)\big ] \end{eqnarray} (B11)with   \begin{eqnarray} R_n = \frac{\Re [K_V(i\, \omega _n)]}{\pi \, n} \end{eqnarray} (B12a)  \begin{eqnarray} I_n = \frac{\Im [K_V(i\, \omega _n)]}{\pi \, n} \end{eqnarray} (B12b)and C(t) is the inverse Laplace transform of $$\tilde{c}(s)/s$$, which can be obtained by integration in the Laplace domain over a closed contour Γ, containing all its singularities,   \begin{eqnarray} C(t) = \mathcal {L}^{-1}[\tilde{c}(s)/s] = \frac{1}{2\, \pi \, i}\, \int _{\Gamma } \frac{\tilde{c}(s)}{s}\, e^{s\, t}\, \mathrm{d}s. \end{eqnarray} (B13)We note that, as depicted in Fig. B1, the closed contour Γ must contain the origin s = 0 and the interval $$(-1/\tau _{\rm {min}},0)\subset \mathbb {R}^-$$ of the negative half of the real axis of the Laplace domain. As implied by the sine and cosine functions in eq. (B11), P(t) is a periodic function with period T and yields zero at t = n T. On the contrary, due to the singularities of $$\tilde{c}(s)$$ having negative real part and zero imaginary part, C(t) has an asymptotic limit for t → ∞, say C∞, that we can estimate by taking the limit of $$\tilde{c}(s)$$ for s → ∞   \begin{eqnarray} C_\infty &=& \lim _{t\rightarrow \infty } C(t) = \lim _{s\rightarrow 0} \tilde{c}(s) = -\frac{\tilde{K}_V(0)}{2} - \tilde{p}(0) \nonumber \\ &=& -\frac{\tilde{K}_V(0)}{2}+\sum _{n=1}^{\infty } I_n \end{eqnarray} (B14)by making use of eqs (B6), (B8) and (B9). APPENDIX C: DISCRETE RELAXATION SPECTRUM The expression obtained in Appendix A greatly simplifies when we consider the case in which the relaxation spectrum of the Maxwell Earth model is discrete. In this case, the viscoelastic Green function takes the following form:   \begin{eqnarray} K(t) = K_E\, \delta (t)+ H(t)\, \sum _{j}k_j\, e^{-t/T_j}\!, \end{eqnarray} (C1)where δ(t) is the Dirac delta, and kj and Tj are the residue and the characteristic relaxation time of the jth relaxation mode. This form, in the Laplace domain, becomes   \begin{eqnarray} \tilde{K}(s) = K_E+\sum _j\frac{k_j}{s-s_j} \end{eqnarray} (C2)with sj = −1/Tj being the pole of the jth relaxation mode, confined on the interval $$[-1/\tau _{\rm {min}},0]\subset \mathbb {R}^-$$ of the negative half of the real axis of the Laplace domain. In this case and considering the dip-slip time history describing the infinite earthquake sequence and back-slip, making use of the residue theorem, the geoid anomaly becomes   \begin{eqnarray} {G(t)= U\, \sum _j\bigg [K_E\, f(t)+r_j\, \left(\frac{1-e^{-t/T_j}}{T/T_j}-\frac{t}{T}\right)} \nonumber \\ {\qquad\qquad\qquad+ \sum _{n=1}^{\infty } r_j\, \big (1-e^{-(t-n\, T)/T_j}\big )\, H(t-n\, T)\bigg ] } \end{eqnarray} (C3)where rj = −kj/sj = kj Tj is the so called strength of the jth relaxation mode, that is the contribution of the mode to the geoid anomaly due to a stick-slip earthquake in the fluid limit, for t → ∞. The first term in the square brackets of the RHS of eq. (C3) is the instantaneous elastic response, the second term is the response to the back-slip due to only stress relaxation and the third term is the summation of the post-seismic response to each earthquake. Considering a time during the nth seismic cycle, t + n T with t ∈ [0, T), and using the expression for finite geometric series for dealing with the post-seismic contribution of the first nth earthquakes, after some algebra, the geoid anomaly during the nth seismic cycle, eq. (10), becomes   \begin{eqnarray} {G_n(t) = U\, \bigg [-K_E\, \frac{t}{T} +\sum _j r_j \bigg (\frac{1-e^{-(t+n\, T)/T_j}}{T/T_j} -\frac{t}{T}} \nonumber \\ {\qquad\qquad\qquad- e^{(T-t)/T_j}\, \frac{1-e^{-n\, T/T_j}}{1-e^{-T/T_j}}\bigg )\bigg ].} \end{eqnarray} (C4)Then, reorganizing this result as in eq. (10), we obtain specific analytical expressions for P(t) and Sn(t), with E(t) given by eq. (11a). Particularly, we have   \begin{eqnarray} P(t) = \sum _jr_j\, p_j(t) \end{eqnarray} (C5a)  \begin{eqnarray} S_n(t) = \sum _j r_j \left(1-e^{-(n\, T+t)/T_j}\right)\left(\frac{T_j}{T}-\frac{1}{1-e^{-T/T_j}}\right) - K_E \end{eqnarray} (C5b)with   \begin{eqnarray} p_j(t) = \frac{1-e^{s_j\, t}}{1-e^{s_j\, T}}-\frac{t}{T}. \end{eqnarray} (C6)We note that   \begin{eqnarray} S_\infty = \sum _j r_j\, \left(\frac{T_j}{T}-\frac{1}{1-e^{-T/T_j}}\right)-K_E. \end{eqnarray} (C7) It can be verified that these expressions are consistent with eqs (B11) and (B14), considering that, in view of eq. (C2), eq. (B12) becomes   \begin{eqnarray} R_n = \frac{1}{\pi \, n}\sum _j\frac{r_j}{1+(T_j\, \omega _n)^2} \end{eqnarray} (C8a)  \begin{eqnarray} I_n = -\sum _j\frac{2\, r_j\, T_j/T}{1+(T_j\, \omega _n^2)} \end{eqnarray} (C8b)and we have   \begin{eqnarray} \tilde{K}_V(0) = \sum _j r_j. \end{eqnarray} (C9) © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations