A sea level equation for seismic perturbations

earthquakes are a potentially important source of relative sea level variations, since they can drive global deformation and simultaneously perturb the gravity ﬁeld of the Earth. For the ﬁrst time, we formalize a gravitationally self-consistent, integral sea level equation suitable for earthquakes, in which we account both for direct effects by the seismic dislocation and for the feedback from water loading associated with sea level changes. Our approach builds upon the well-established theory ﬁrst proposed in the realm of glacio-isostatic adjustment modelling. The seismic sea level equation is numerically implemented to model sea level signals following the 2004 Sumatra–Andaman earthquake, showing that surface loading from ocean water redistribution (so far ignored in post-seismic deformation modelling) may account for a signiﬁcant fraction of the total computed post-seismic sea level variation.


I N T RO D U C T I O N
The problem of future climate changes and their impact on human activities is still far from a complete solution (IPCC 2007). Nevertheless, the computational efforts devoted to this problem have reached a phase in which second-order complexities are often taken into account to achieve real world resolution levels (Shukla et al. 2006;Bamber et al. 2009;Mitrovica et al. 2009). The sea level variation associated with seismic activity is a representative of these second-order effects. In fact, long term sea level change is driven primarily by eustatic processes, glacio-isostatic adjustments (GIAs) and thermocompositional volume variations (IPCC 2007), while seismo-tectonic deformations play a small (but non negligible) role (Melini et al. 2004;Melini & Piersanti 2006).
The main shortcoming of the investigations so far is that, instead of approaching the full sea level equation (SLE), an approximated solution was computed only taking into account the direct effects of seismic sources on deformation and gravity potential variations (e.g. Melini & Piersanti 2006). This is justified in regional postseismic investigations on timescales of a few decades for which the isostatic response is expected to be negligible, but it is likely to be incorrect in global studies, where self-gravitation of the oceans plays a fundamental role (Farrell & Clark 1976). Recently, De Linage et al. (2009) solved a zeroth-order sea level equation for the short-term relaxation following the 2004 Sumatra-Andaman earthquake; according to their results, the response of the ocean has to be taken into account in order to correctly interpret the observed geoid perturbations. To improve these aspects of post-seismic rebound modelling, in this short note we describe and numerically solve a gravitationally self-consistent SLE for seismic perturbations, generalizing the results of Farrell & Clark (1976). In our study, the post-seismic deformation and gravity potential variation are first obtained by a semi-analytical, self-gravitating viscoelastic model (Piersanti et al. 1995), and are then used as initial conditions for an iterative solution scheme for the SLE, in which the loading problem associated with the mass redistribution of the oceans is solved using a post-glacial rebound calculator .
In Section 2, we discuss the theoretical aspects of our method and in Section 3 we focus on the details of numerical implementation, investigating the convergence of the iterative solution and addressing a simple synthetic problem. In Section 4, the method is used to evaluate the global and regional sea level variations following the Sumatra-Andaman earthquake of 2004 December 26.

M E T H O D S
According to the theory of Farrell & Clark (1976), in the framework of GIA the SLE reads where S is sea level change, ω = (θ, λ) denotes colatitude and longitude, t is time, is the incremental gravity potential, γ is the reference gravity acceleration at the Earth surface, U is vertical displacement and S E is the eustatic sea level change, which represents the solution of the SLE in the case of a rigid, non-self-gravitating Earth where m i is the mass variation of ice sheets, ρ w is the density of water, A o = oceans d A denotes the (constant) area of the ocean Sea level equation for seismic perturbations 89 surface and the overline indicates the average A SLE suitable for seismic perturbations can be readily obtained from eq. (1) dropping the S E term, since earthquakes do not imply any mass exchange with the oceans. However, the averaged term must be kept, since it ensures that the free surface of the oceans always coincides with the geoid (Farrell & Clark 1976). This gives where S(ω, t) defines the history of sea level change at any point ω = (θ, λ) on the sphere and where now and U are the total gravity potential variation and surface displacement imposed by the seismic dislocation, respectively. Consistently with the principle of mass conservation, in the seismic SLE (4), S = 0.
The terms U and in eq. (4) stem from the sum of two contributions. The first, labelled by 'eq' below, represents the direct effect of the seismic dislocation, while the second ('load') is associated with the water load exerted by the changing sea level. Such decomposition is similar to that adopted in the framework of GIA studies (Spada & Stocchi 2006). Thus we write  Tegmark (1996) pixelization. The dots are the centroids of slightly distorted, equal-area hexagons. The code that generates the pixels coordinates is available from http://space.mit.edu/home/tegmark/. Wet and dry pixels are separated from the global distribution using the GMT utilities (Wessel & Smith 1991), freely downloadable from http://gmt.soest.hawaii.edu/. 90 D. Melini, G. Spada and A. Piersanti and (ω, t) where the S-dependence of the load terms can be expressed by a time-convolution between the viscoelastic loading-deformation coefficients h l (t) and k l (t) and the history of sea level change, which makes eq. (4) an integral (implicit) equation. The lack of the eustatic term and the simple Heaviside time-history usually employed to describe the source (e.g. Piersanti et al. 1995) makes the seismic SLE formally simpler but does not alleviate the numerical complexity of the problem, as will be discussed below. In previous studies (Melini et al. 2004;Melini & Piersanti 2006), the ocean-averaged term in eq. (4) was neglected. For an incompressible Earth, this is equivalent to the assumption of a uniform, non-self-gravitating ocean. Furthermore, the approximation load = U load = 0 was adopted, which reduces the SLE to a fully explicit equation that can be solved as soon as the direct effect of earthquakes is determined. A zeroth-order approximation to the solution S of the SLE (4) can be obtained neglecting load and U load in front of the 'eq' terms. With eq ≡ (0) and U eq ≡ U (0) , this gives which is used to provide a first guess of the water load (mass per unit area) according to where O is the 'ocean function' (O = 1 over the oceans, and O = 0 elsewhere) and where positive and negative values of L correspond to a sea level rise and fall, respectively. Once L (0) is determined globally, the response to loading U load and load can be computed using pertinent load-deformation coefficients, providing a new estimate of the total displacement and gravity potential (1) (ω, t) = eq + load (S (0) ), which substituted into the right-hand side of eq. (4) gives the new estimate of sea level change, S (1) . The method outlined above suggests the following general algorithm: (i) given S (k) , the kth order approximation of the sea level change (k = 0, 1, . . .), compute the water load function L (k) (ω, t) = ρ w S (k) O by eq. (9), (ii) using the direct responses to seismic dislocation and the solution to the loading problem, evaluate U (k+1) = U eq + U load (S (k) ) and (k+1) = eq + load (S (k) ), (iii) from the SLE (4), compute the further approximation to sea level change S (k+1) , (iv) iterate until a previously defined convergence criterion is satisfied and (v) if needed, provide final estimates for the total perturbation to gravity potential and vertical displacement field. This scheme is largely similar to that employed in GIA investigations, which has been thoroughly validated in a number of case-studies (see Spada & Stocchi 2006, and references therein), generally showing a fast convergence.

N U M E R I C A L I M P L E M E N TAT I O N
In our implementation, the response functions U eq and eq in eqs (5) and (6) are computed by the viscoelastic normal-mode approach originally proposed by Piersanti et al. (1995), for an incompressible, spherical self-gravitating model with Maxwell rheology. The algorithm outlined in Section 2 could be also applied to finely layered earth models, possibly characterized by a generalized (linear) rheology (Spada & Boschi 2006;Melini et al. 2008;Spada 2008) or mantle compressibility (Pollitz 1997(Pollitz , 2003. The response of the Earth to surface loading is evaluated by suitably adapting the TABOO post-glacial rebound calculator , see http://samizdat.mines.edu/taboo/). The model, described in Table 1, is characterized by a coarse 4-layer structure with PREM-averaged density and rigidity and includes a low-viscosity upper mantle beneath a perfectly elastic lithosphere, and an homogeneous inviscid core. The relatively low value of upper mantle viscosity has been chosen because of the importance of the low-viscosity zone in postseismic rebound processes (Nostro et al. 1999;Piersanti et al. 2001). Since the solutions obtained by Piersanti et al. (1995) are limited to a four-layer stratification, we decided to employ for the upper mantle a viscosity that could reproduce effects similar to those associated with the presence of a low-viscosity zone.
The computation of surface integrals in eq. (4) and of the responses to water load (i.e. functions U load and load in eqs (5) and (6)) are practically performed with the aid of the icosahedron-based pixelization shown in Fig. 1, proposed by Tegmark (1996) for astrophysical applications. In the GIA context, this grid has been employed for the first time by Spada & Stocchi (2007) for solving numerically the SLE. The Tegmark discretization provides a natural set of Gauss points on the surface of the sphere and allows for a straightforward computation of surface integrals involving spherical harmonic (SH) functions as equal-weight finite sums. This property can be employed to compute ocean-averages as where f is a scalar function, N is the total number of pixels (according to Tegmark 1996, N = 40R(R − 1) + 12 where R is a resolution parameter), N w is the number of ocean pixels, ω i are their coordinates and O 00 0.71 is the degree zero and order zero harmonic coefficient of the ocean function (4π -normalized complex SH will be used throughout). For a given grid resolution R, Tegmark (1996) has shown that approximation (12) is numerically valid as long as the maximum degree of the SH expansion of f is l max ≤ √ 3N . The icosahedron-based pixelization is also employed to discretize the surface load defined by eq. (9). The load is distributed over axisymmetrical disc-shaped elements with centres defined by the ocean pixels of Fig. 1, each with a diameter d = 2a arccos(1 − 2 N ), a being Earth radius. Since resolving each load component requires an SH expansion to a degree l max 2πa/d, a correct numerical implementation of the SLE thus requires π   Table 2). that allows an optimal trade-off between grid spacing and computational costs to be determined. To satisfy eq. (13), in our simulations we have used a grid with N = 15212 (this corresponds to R = 20) and considered harmonic degrees up to l max = 200. The computation of U load and load takes advantage of the symmetry of the load components, which makes these terms only dependent upon the colatitude of the observer relative to each elementary disc. The convolution integrals that involve load-deformation coefficients h l (t) and k l (t) and the history of sea level change within each disc load are discretized in the time domain and computed by standard numerical methods .
To test the stability and the convergence of the solution scheme and to verify the absence of aliasing effects due to pixelization, we have performed a test imposing ad hoc seismic effects. In particular, we set eq = 0 and prescribe, for time t ≥ 0, a vertical displacement U eq = −1 m across a circular region of half-amplitude α = 5 • placed at ω = (π/2, π ) (i.e. in the central Pacific Ocean). Thus, in this experiment, where U eq is negligible in front of U eq because of the localized displacement assumed for our toy earthquake, and using eq. (4) the load-induced sea level variation is In Fig. 2, S load is shown for k = 1, 2, 5 and 10, as a function of time and for various source-observer angular distances . For = 0 • and relatively short times, S (k) load slightly enhances the direct seismic effect S (k) eq ≡ −U eq . However, with increasing time, S load becomes a large fraction (∼40 per cent) of the direct effect in the vicinity of the seismically deformed region, due to the viscoelastic relaxation induced by the water load. S (k) load is large in the vicinity of the source and decays quickly with the observer distance, falling by a factor of ∼10 2 moving from = 0 (frame a) to = 20 • (d). It is interesting to observe that, in spite of the low-viscosity upper mantle (see Table 1), the S (k) load curves are still far from equilibrium at time t = 1 kyr after loading, a timescale that largely exceeds the Maxwell relaxation time for this layer (3.7 yr). This may be interpreted as an effect of the response of the lower mantle, which is involved due to the relatively large area of the 'fault plane' employed in this synthetic case study. The density jump imposed at the depth of 670 km is also likely to play a role, due to the long relaxation times that characterize the return of compositional boundaries to equilibrium (Piersanti et al. 1995).
The issue of the convergence of the iterative scheme is addressed more quantitatively in Fig. 3, where the ratio S load is shown as a function of k for = 0 (a) and = 20 • (b), and various times following the synthetic earthquake already considered in the previous figure. It is apparent that the convergence is monotonic and relatively fast: these features are qualitatively similar to those observed when the SLE is solved for glacial forcing (e.g. Spada & Stocchi 2007). The spherically averaged relative difference between subsequent iterations, defined as (k)  Fig. 4 for various values of time t, indicates that k = 5 ensures incremental errors well below the 0.1 per cent threshold, which is fully acceptable for any practical implementation.

S E A L E V E L VA R I AT I O N S F O L L O W I N G T H E 2 0 0 S U M AT R A -A N DA M A N E A RT H Q UA K E
In this section, we present an application of the proposed method to the sea level variations following the 2004 Sumatra-Andaman earthquake. The seismic source has been modelled with five point dislocations corresponding to the multiple CMT solution obtained by Tsai et al. (2005). These sources are obtained by fitting with the CMT method the long-period seismograms from the IRIS Global Seismographic Network. They account for a cumulative energy release corresponding to M w = 9.3; their locations and focal mechanisms are shown in Fig. 5. Using the semi-analytical model of global post-seismic rebound originally developed by Piersanti et al. (1995) and subsequently extended by Soldati et al. (1998) and Boschi et al. (2000), we have obtained the time-dependent deformation and incremental gravitational potential U eq and eq . These fields have been used as starting conditions to iteratively solve eq. (4), as discussed in Section 2.
To evaluate the zeroth-order solution of the SLE defined in eq. (8), we need to compute oceanic averages of U eq and eq according to eq. (3). Since the body-force equivalent representation of a point source is based on localized Dirac delta functions and their spatial derivatives (Smylie & Mansinha 1971;Mansinha et al. 1979), the spectra of U eq and eq are rich at short wavelengths and decay slowly with harmonic degree (Casarotti 2003;Melini et al. 2008). For this reason, to obtain convergence, the SH expansion of the relevant scalar fields (U eq and eq ) has to be truncated to l max ≈ 10 3 -10 4 , depending on the source-observer distance (Sun & Okubo 1993;Riva & Vermeersen 2002;Casarotti 2003). In the present application, the post-seismic solutions reach a stable convergence for l max = 4000. Since for numerical stability of eq. (3) the relation l max ≤ √ 3N must be satisfied, the computation of U eq and eq requires a pixelization with N ≥ 5.3 × 10 6 points, corresponding to a resolution R = 366. We remark that, due to the linearity of oceanic integrals, this high-resolution pixelization is not needed in subsequent iterations of the SLE solution scheme. Indeed, once U eq and eq are known, the evaluation of oceanic integrals of eqs (5) and (6) requires only the integration of 'load' terms, which can safely be carried out with an R = 20 pixelization, as discussed in the previous section. The computation of U eq and eq on the high-resolution grid represents a very intensive numerical task: even with a highly optimized parallel integration code on a 128core distributed-memory cluster, it requires about 15 hr for each point source. For this reason, while a 2-D source modelling would be certainly more realistic (Nostro et al. 1999), we are currently limited to the point-source approximation; indeed, considering the exceptionally large rupture extension of the Sumatra earthquake, modelling a 2-D source through a superposition of point sources would increase the computation time above acceptable levels, even adopting a relatively coarse source discretization. We remark that the iterative solution method for the SLE, which has been outlined in Section 2, is independent from the post-seismic deformation model, which is used only to provide initial conditions to the iterative solution scheme. The seismic SLE can therefore be solved with the same prescriptions if a more detailed model of post-seismic deformation is employed. For the present application, whose aim is a Figure 8. As in Fig. 7, but on a regional scale. demonstration of the water load effects in a real case, we will therefore use a point-source approximation, which anyway gives acceptable results on a global scale.
In Fig. 6, the average relative difference between iterations (k) r (t), as defined by eq. (16), is shown for a range of observation times. From a comparison with Fig. 4 it can be observed that, when a real seismic source is employed, the convergence of the iteration scheme is less regular than in the synthetic case. This is likely to be the result of the increased numerical noise introduced by the rich spectrum of harmonics that characterize the realistic seismic source compared with the 'hat' test displacement considered in Section 3. In spite of this, however, after k = 4 iterations the average relative difference is ≤ 5 per cent and for k = 10 it is below the 1 per cent level. Looking at the spatial patterns of S (k) (ω, t), we have verified that less regular convergence specifically results from contributions to (k) r (t) from regions close to the nodal lines of this function, where some of the terms in eq. (16) become numerically indeterminate, because of S (k) 0, even if the solution has already reached a stable convergence in the bulk of the spatial domain.
In Figs 7 and 8, we quantitatively evaluate the effect of water load upon sea level changes, focussing on a global and a regional scale, respectively. The left-hand frames show snapshots of S (k=0) (ω, t), which only accounts for the effect of the seismic dislocation source, computed according to eq. (8), while in the right-hand frames we consider S load = S (k=10) − S (k=0) . From Fig. 7, the term S load turns out to be smaller than the seismic contribution, but definitely not negligible, being a significant fraction of the total signal even on a global scale. Its relative weight increases with time, due to the delayed viscoelastic response of the ductile layers to the forcing of the seismically induced sea level variations. For the local scale analysis of Fig. 8, it results that the effect of the water load correction is even stronger, but the results may be biased to some degree by the point-source approximation which can affect significantly the near-field computations (Nostro et al. 1999). For short timescales (a few years) the load correction is manifest as a broad sea level fall, with a smoothed pattern with respect to the negative lobe associated with the purely seismic contribution, consistently with the results of the synthetic case discussed above. For longer timescales (t = 100 yr in Fig. 8) the contribution to sea level from water load broadly follows the pattern of alternating lobes of the seismic term; this is a common feature found in post-seismic relaxation of low-dip thrust faults (Rundle 1982;Volpe et al. 2007).  Figure 10. The same as in Fig. 9, but on a time window of 20 years following the main shock. Dotted lines show least-squares trends of the individual solutions. Fig. 9 shows predictions of post-seismic sea level variations at ten tide-gauge sites belonging to the PSMSL network (see http://www.pol.ac.uk/psmsl/), whose locations are shown in Fig. 5. Solid and dashed curves show results obtained including and neglecting the water load in the SLE, respectively. All the sites share a qualitatively similar history of sea level variations, in which a postseismic viscoelastic 'wave' follows a quiescent initial phase on timescales of a few centuries. As we have found with the synthetic tests of Fig. 2, this timescale exceeds the intrinsic Maxwell characteristic time of mantle layers, suggesting that relaxation of internal compositional boundaries are indeed playing a role, due to the large extent of the seismic source (Piersanti et al. 1995). At near-field sites (e.g. Ko Taphao Noi), our computations predict a sea level fall of ∼1 m during the next century, which would imply average rates of sea level change that greatly exceed the secular globally averaged trend, close to 1.5 mm yr −1 (IPCC 2007). Sensibly smaller (but still significant) effects are predicted for other sites (Kanmen, Manila and Danang), with a sea level rise of up to ∼5 cm during the same period. The average trend, in this case, is ∼30 per cent of the current average global trend, and practically negligible in comparison with the local sea level trend which amounts to ∼12 mm yr −1 during the last few decades in the case of Manila (Spencer & Woodworth 1993). While corrected and uncorrected sea level predictions are generally similar on a decade timescale, they may diverge for longer periods, when the water loading effect may perturb the seismic contribution significantly, by values ranging between 10 and 20 per cent. This is also found for 'far-field' tide-gauges (e.g. Port Louis and Broome), where numerical artefacts due to the point-source approximation are likely to have a minor role. In Fig. 10 the synthetic sea level time-series are plotted on a 20-yr period, during which they are well approximated by a linear trend. With a least-squares linear regression (see dotted lines), an estimate of the rate of sea level variation has been obtained from the results of Fig. 10; numerical values are listed in Table 2. The contribution to the rate of sea level variation due to the load correction,Ṡ load , turns out to be a large fraction the total trend, with a relative impact up to nearly 50 per cent at Vishakhapatnam. The coseismic and post-seismic gravity field perturbations following the 2004 Sumatra earthquake have been evidenced by GRACE satellite measurements. Several authors have extracted the earthquake signature from the GRACE solutions and found it to be consistent with seismological models (Han et al. 2006;Ogawa & Heki 2007;Panet et al. 2007). Recently, De Linage et al. (2009) modelled the post-seismic geoid perturbation taking into account the static potential perturbation of a global incompressible ocean; according to their results, the oceanic contribution is needed in order to successfully reproduce the spatial features of the post-seismic geoid perturbation observed by GRACE. As we verified by extracting the geoid signal from our results, the oceanic term obtained by De Linage et al. (2009) has the same sign and spatial extension as the short-term water load effect resulting from our simulations, even if the seismic source model employed by De Linage et al.
(2009) is more realistic than our point-source approximation. Our simulations predict larger peak values of geoid perturbation, which is probably a bias effect of the point source model which leads to an overestimation of coseismic effects in the near-field (Piersanti et al. 1997).
The spectral features of the sea level correction due to ocean loading can be investigated by computing the harmonic coefficients Figure 11. Non-dimensional, normalized spectral coefficients in the range of harmonic degrees 2 ≤ l ≤ 120 computed according to eq. (17).
where is the unit sphere and Y lm are the 4π-normalized complex SH. The normalized squared coefficients |c lm | 2 / max l,m (|c lm | 2 ), displayed in the diagrams of Fig. 11, show that on timescales of a few years most of the signal is confined to low harmonic degrees, while for longer times the relative weight of higher harmonics increases. This indicates the presence of short-wavelength features of the sea level signal in the area surrounding the seismic source as a consequence of stress concentration due to viscoelastic relaxation in the ductile layer, which for low-angle thrust faults may result in smallscale regions of opposite vertical deformation around the source location (Rundle 1982;Volpe et al. 2007).

C O N C L U S I O N S
For the first time we have obtained a solution for the gravitationally self-consistent SLE describing sea level perturbations occurring after a large earthquake. The SLE has been solved numerically by implementing an iterative scheme directly derived from those adopted in GIA studies (see e.g. Spada & Stocchi 2006). As a result, our analysis shows that feedback loading effects play a significant role in assessing seismic quasi-static sea level variations. The viability of the proposed approach has been assessed by means of a synthetic test with a disc-shaped oceanic load in order to verify its numerical stability. The solution convergence turns out to be monotonic and relatively fast, similarly to what is observed in postglacial rebound applications (Spada & Stocchi 2007). The method has then been applied to the prediction of sea level variations following the 2004 Sumatra-Andaman earthquake. We found that loading effects represent an important contribution to seismically induced sea level variations on timescales ranging from a few decades up to several thousands of years. These timescales, which largely exceed the Maxwell relaxation times of the involved layers, suggest that relaxation modes connected to internal compositional boundaries are excited. An analysis of the predicted sea level signal on a set of PSMSL tide-gauge sites showed that, for 'nearfield' stations, the expected post-seismic effect is not negligible even in comparison with the globally averaged secular trend, although this result may be biased by the point-source approximation which is currently unavoidable due to computational requirements.
The presence of long-term effects suggests that a detailed knowledge of historical seismicity is crucial in modelling present-day sea level rates. For timescales of a few years, the sea level signal follows an approximately linear trend, and the loading term represents a non-negligible perturbation to the total rate. These short-term effects may be further enhanced in the presence of rheological layers characterized by a transient rheology, since in that case a large postseismic signal occurs on timescales of the order of months (Pollitz 2003).
In this respect, we can conclude that a detailed modelling of sea level change cannot neglect the effect of seismic perturbations, which can be the predominant contribution in correspondence of subduction zones characterized by large seismic energy release. Future high-resolution scenarios of sea level variation should take into account, among other contributions, the highly heterogeneous signals coming from short wavelength regional seismic activity, in order to precisely assess the exact role played by different phenomena in determining sea level variation. The inclusion of seismic effects in a comprehensive approach based on a self-consistent solution of the SLE represents an opportunity to create a unified formal framework to model non-eustatic sea level variations.

A C K N O W L E D G M E N T S
We thank John Beavan, Wenke Sun and Riccardo Riva for their careful and incisive reviews. This work is partly supported by MIUR (Ministero dell'Istruzione, dell'Università e della Ricerca) with the PRIN2006 grant 'Il ruolo del riaggiustamento isostatico postglaciale nelle variazioni del livello marino globale e mediterraneo: nuovi vincoli geofisici, geologici ed archeologici' and the FIRB grant 'Sviluppo di nuove tecnologie per la protezione e la difesa del territorio dai rischi naturali'. The figures have been prepared with the GMT software by Wessel & Smith (1991). The post glacial rebound solver TABOO is available from the Samizdat Press (http://samizdat.mines.edu/taboo/).