Portal de Periódicos da CAPES

A model for seasonal changes in GPS positions and seismic wave speeds due to thermoelastic and hydrologic variations

2011; American Geophysical Union; Volume: 116; Issue: B4 Linguagem: Inglês

10.1029/2010jb008156

2156-2202

Victor C. Tsai,

High-pressure geophysics and materials

Journal of Geophysical Research: Solid EarthVolume 116, Issue B4 Geodesy and Gravity/TectonophysicsFree Access A model for seasonal changes in GPS positions and seismic wave speeds due to thermoelastic and hydrologic variations Victor C. Tsai, Victor C. Tsai [email protected] Geologic Hazards Science Center, U.S. Geological Survey, Golden, Colorado, USASearch for more papers by this author Victor C. Tsai, Victor C. Tsai [email protected] Geologic Hazards Science Center, U.S. Geological Survey, Golden, Colorado, USASearch for more papers by this author First published: 19 April 2011 https://doi.org/10.1029/2010JB008156Citations: 89AboutSectionsPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Abstract [1] It is known that GPS time series contain a seasonal variation that is not due to tectonic motions, and it has recently been shown that crustal seismic velocities may also vary seasonally. In order to explain these changes, a number of hypotheses have been given, among which thermoelastic and hydrology-induced stresses and strains are leading candidates. Unfortunately, though, since a general framework does not exist for understanding such seasonal variations, it is currently not possible to quickly evaluate the plausibility of these hypotheses. To fill this gap in the literature, I generalize a two-dimensional thermoelastic strain model to provide an analytic solution for the displacements and wave speed changes due to either thermoelastic stresses or hydrologic loading, which consists of poroelastic stresses and purely elastic stresses. The thermoelastic model assumes a periodic surface temperature, and the hydrologic models similarly assume a periodic near-surface water load. Since all three models are two-dimensional and periodic, they are expected to only approximate any realistic scenario; but the models nonetheless provide a quantitative framework for estimating the effects of thermoelastic and hydrologic variations. Quantitative comparison between the models and observations is further complicated by the large uncertainty in some of the relevant parameters. Despite this uncertainty, though, I find that maximum realistic thermoelastic effects are unlikely to explain a large fraction of the observed annual variation in a typical GPS displacement time series or of the observed annual variations in seismic wave speeds in southern California. Hydrologic loading, on the other hand, may be able to explain a larger fraction of both the annual variations in displacements and seismic wave speeds. Neither model is likely to explain all of the seismic wave speed variations inferred from observations. However, more definitive conclusions cannot be made until the model parameters are better constrained. Key Points To model seasonal variability in GPS displacements and seismic wave speeds To model this variability with a thermoelastic model and a hydrologic model To tentatively compare the model results with observations 1. Introduction [2] Due to the importance of thermoelastic stresses and strains in engineered materials, the theory of thermoelasticity has a long and well developed history [see, e.g., Love, 1944; Timoshenko and Goodier, 1970]. Only in the past few decades, however, have thermoelastic effects been examined in a geologic context. In particular, Berger [1975] developed the theory to describe thermoelastic strains in a two-dimensional (2-D) homogeneous elastic half-space with a periodic surface temperature variation, and used this theory to describe how observed periodic strains within the Earth can be modeled satisfactorily. Ben-Zion and Leary [1986] expanded upon Berger's work to explicitly treat the same case but with an unconsolidated, incompetent surface layer above the elastic half-space. They found reasonable agreement between observed strains at a few sites in southern California and strains modeled using observed temperature variations, with typical strains being on the order of 10−7–10−6 and being delayed relative to the temperature variation by 2 to 3 months. [3] Hydrologic variability has also been studied for many years [e.g., Todd, 1959], primarily due to the importance of groundwater as a source of drinking water. However, the loading on the Earth produced by these variations has not been studied until more recently [e.g., Crowley et al., 2006; Bettinelli et al., 2008]. Fortunately, the theory that governs both the direct elastic loading and the poroelastic loading due to hydrologic variability is well understood. In fact, the theory of poroelasticity is known to have a one-to-one correspondence with the theory of thermoelasticity [Biot, 1956; Rice and Cleary, 1976]. [4] GPS time series are well known to be affected by a potentially large number of factors, including tectonic motions, tidal loading, and the hydrologic cycle [Dong et al., 2002]. Prawirodirdjo et al. [2006] further suggested that thermoelastic strains contribute significantly to the observed seasonal periodicity in GPS displacements (which has typical annual variations on the order of a few millimeters). Prawirodirdjo et al. [2006] found that the phase of the yearly variations in GPS positions could be well matched using the thermoelastic model of Ben-Zion and Leary [1986] with reasonable choices of parameters. However, in this study, the amplitudes of these variations were not considered, even though the displacements expected from thermoelastic strains are easily calculated by integrating the strains. In section 2.4, I provide expressions for these displacements and compare the calculated amplitudes with those observed (section 3.1). I find that estimates of the maximum realistic thermoelastic displacements can account for a fraction of the observed GPS signal at a typical site in southern California. Bettinelli et al. [2008] have also shown that direct elastic loading from water table variations can contribute to significant GPS displacements, and in section 2.4 I also calculate displacements corresponding with both this elastic effect and the poroelastic effect. I find that the maximum realistic hydrology-induced displacements can also account for a significant fraction of the observed GPS signal. [5] Even more recently, by examining the coda of cross correlations of ambient seismic noise, Meier et al. [2010] have found a seasonal variation in observed seismic travel times in southern California. While there remains some possibility that these observations are due to seasonal variations in the location of seismic noise sources [e.g., Zhan and Clayton, 2010], the most straightforward interpretation of the observations is in terms of seasonal changes in seismic velocity structure on the order of δV/V ≈ 0.1%. Meier et al. [2010] suggest that two possible reasons for these temporal variations in seismic wave speed are hydrologic and/or thermoelastic variations, but do not estimate whether the expected variability can produce the observed results. In section 2.5, I evaluate both of these hypotheses by calculating the expected wave speed variations using the third-order elasticity theory of Murnaghan [1951]. After performing order-of-magnitude estimates of certain parameters, I find that the calculated wave speed variations are unlikely to explain most of the observed variations (section 3.2), but that uncertainties in the parameters precludes making a stronger conclusion. [6] Finally, before turning to the details, I would like to stress that the model results presented here, like many purely analytical results, can only be expected to approximate any given realistic scenario. Not only will it be seen that some of the model parameters are poorly constrained, but a number of approximations are also used that are not expected to hold perfectly. The usefulness of these analytical results are in understanding how quantitative results depend on the various key physical parameters, and in allowing researchers to quickly evaluate the plausibility of certain models explaining the observed phenomena. Since such plausibility studies have not yet been done to the extent performed here, I believe it to be useful for the community, as a step toward understanding the observations. I further believe that attempting more accurate and/or more precise quantitative analysis is not currently prudent, and should await further constraints. 2. Theoretical Results in Thermoelasticity and Poroelasticity Berger's Thermoelastic Strain Solution [7] Berger [1975] derived a solution to the uncoupled quasistatic thermoelastic equations [e.g., Timoshenko and Goodier, 1970; Boley and Weiner, 1997] for a 2-D elastic half-space in plane strain subjected to a sinusoidal surface temperature variation in both space and time. This 2-D approximation to a truly 3-D situation is used throughout this work; as in other periodic elasticity problems, 3-D effects are of second-order importance especially when wavelengths in one direction are considerably shorter than in the other direction [Malvern, 1969], and I do not attempt to include any of these dependencies. With x denoting horizontal position, y denoting depth, and t denoting time, the surface temperature boundary condition is given by Tavg is a constant background average temperature, T0 is half of the peak-to-peak amplitude of the periodic temperature variation, k is the horizontal wave number, and ω is frequency (see Figure 1). Assuming the temperature variation approaches zero at depth (T(x, y → ∞, t) − Tavg = 0), Berger finds expressions for the strains εxx, εyy and εxy. These expressions contain "equivalent body force" terms that decay with depth in proportion to the temperature variation and "surface traction" terms that decay with depth proportional to e−ky. Figure 1Open in figure viewerPowerPoint Schematic of the thermoelastic model. TS is the surface temperature variation, yT is the thermal boundary layer thickness, yb is the thickness of the incompetent layer, and L is the wavelength of the surface temperature variation. The GPS station (with position ux and uy) is assumed anchored below yT, and seismic velocities, Vij, are calculated in the elastic half-space. Schematic depth profiles for temperature, T(y), and strain, εij(y), are graphed. [8] For annual variations on the Earth (ω ≈ 2 · 10−7 s−1) with kilometer-length scales and thermal diffusivity ≈ 10−6 m2/s [Berger, 1975], then κk2/ω ≪ 1, and Berger finds that the expressions simplify significantly, especially for depths not within the thermal boundary layer. Using these approximations and accounting for an incompetent layer of thickness yb, as in the work of Ben-Zion and Leary [1986] (see Figure 1), I obtain simple expressions for the strains and where ν is Poisson's ratio and A(t) is given by In this expression, αth is the coefficient of linear thermal expansion (αth ≈ 10−5°C−1) [Berger, 1975], and equations (2) are appropriate for depths outside the thermal boundary layer, i.e., y ≫ . Examining equations (2), one observes that all components of the thermoelastic strain have the same time dependence, similar amplitudes, and decay with depth in a similar fashion. Equation (3) shows that the strains are delayed relative to the temperature variation by an amount where τ ≡ 2π/ω is the period. [9] Using equation (4) to solve for yb in terms of Δt (and other parameters), we can rewrite A(t) as a function of Δt as The strain amplitude versus phase lag of equation (5) is plotted in Figure 2 with ν = 0.3, k = 2π/(10 km), T0 = 10°C, and other values as before (ω ≈ 2 · 10−7 s−1, κ ≈ 10−6 m2/s). The value of T0 chosen is on the high end of representative temperature variations [Berger, 1975]. Figure 2Open in figure viewerPowerPoint Strain amplitude A(t) (amplitude only, without time dependence) as a function of time delay Δt. Thermoelastic results are for varying yb, poroelastic results are for varying κhy, and water table elastic results are for varying ϕ. Maximum strain and minimum time delay for the thermoelastic solution occur when the incompetent layer thickness, yb, is small. Maximum strain and minimum time delay for the poroelastic solution occur for the largest values of κhy. Maximum strain for the water table elastic solution occurs at the largest values of ϕ. Representative values of parameters are chosen as described in the text (e.g., T0 = 10°C, p0 = 2.9 · 104 Pa). For Figure 2, k = 2π/(10 km). Note that the thermoelastic and poroelastic amplitudes scale directly with k and that in section 3 it is suggested that k may actually be closer to 2π/(40 km). Poroelastic Hydrologic Modification [10] Since thermoelasticity and poroelasticity share the same mathematical framework [Biot, 1956; Rice and Cleary, 1976; Wang, 2000], the corresponding decoupled quasistatic poroelastic problem is solved with the correspondence where E is Young's modulus, α is the Biot-Willis coefficient, p is pore pressure, and κhy is hydraulic diffusivity. For this problem, the boundary condition is where pavg is a constant pressure and p0 is the amplitude of the pore pressure variation. In reality, such a boundary condition does not exist but can be used to approximate variations in water table level. In addition, the decoupled approximation is not completely valid [Detournay and Cheng, 1993], but gives a reasonable approximation to the fully coupled problem [Roeloffs, 1988]. [11] For a nominal hy ≈ 0.4 m2/s, hy ≡ hyk2/ω ≈ 0.8 does not satisfy hy ≪ 1 as in the thermoelastic solution. However, the approximate solution is still roughly valid as long as hy ≲ 1, with modification to equation (4) so that where the approximation is valid for κhy ≳ 10−3 m2/s (and yb < 10 m). (Note that Berger's full unapproximated solution could be used if better precision were desired.) The approximate poroelastic solution is then given by equations (2) and (3) with equation (8) substituted for equation (4) and equations (6) replacing thermoelastic parameters. The strain amplitude versus phase lag is plotted in Figure 2 for a range of κhy, E = 1.6 · 1010 Pa, p0 = 2.9 · 104 Pa, α = 0.7, and other values as before. It is of interest to note that the strain amplitude and phase lag for the poroelastic solution for κhy ≈ 0.05 m2/s is nearly the same as that of the thermoelastic solution for yb = 0. Direct Elastic Hydrologic Modification [12] In addition to the poroelastic effect just discussed, water table fluctuations also produce a direct elastic loading due to the weight of the additional fluid. This purely elastic loading can be calculated with the same methods as given by Berger [1975], yielding where Ae(t) is given by and is porosity (so that fluid weight per unit area is ϕp). The primary difference between these strains and the poroelastic strains is that these strain amplitudes scale with ϕp0/E instead of αp0/E and the strains are exactly in phase with the water level forcing (Δt = 0). The strain amplitude for this direct elastic strain from water table fluctuations is plotted in Figure 2 for a range of ϕ, and other values as before. Note that for ϕ = 0.2, the strain amplitude is similar to the poroelastic strain amplitude when κhy = 0.15 m2/s. GPS Displacements From Modeled Strain [13] Given the strains of equations (2), it is straightforward to integrate to obtain displacements. By assumption, horizontal and vertical displacements, ux and uy, approach zero as y → ∞ and rigid body motion is ignored. Therefore and [14] For a GPS station anchored at a point beneath the thermal boundary layer (see Figure 1), ux and uy of equations (11a) and (11b) are the temporal variations in GPS position expected of thermoelastic variations. GPS stations anchored within the thermal boundary layer (including at the surface) will include additional terms related to terms ignored in the approximations of equations (2). Under the approximations made, these additional terms are only significant for the vertical displacement, uy. Integrating Berger's full solution, one finds that the magnitude of uy would be increased by about a factor of 3 at the near surface if these terms were included. For the purposes of this paper, GPS stations will be assumed to be deeply anchored. [15] For water table fluctuations, the same analysis applies. Equations (11a) and (11b) apply to poroelastic displacements with the modifications of equations (6) and (8), and equations (9) can similarly be integrated to obtain displacements for the purely elastic component. However, since actual water table fluctuations occur at depths of a few meters and furthermore κhy ≫ κ, the GPS stations are likely anchored within the "porous" boundary layer. Thus, there are additional poroelastic displacements as discussed in the previous paragraph. Since only the vertical displacement, uy, is significantly affected and only horizontal displacements are compared, these additional displacements will not be discussed further. Seismic Wave Speed Variation From Modeled Strain [16] For a given strain field, expressions for elastic wave speeds can be obtained by keeping up to third-order terms in the strain energy function [Murnaghan, 1951; Norris, 1998]. For an initially isotropic body, Hughes and Kelly [1953] and Egle and Bray [1976] provide V11, V12 and V13 as and where Vij is the speed propagating in the i direction with polarization in the j direction, ρ0 is the initial density, λ and μ are the Lame constants (λ = Eν/[(1 + ν)(1 − 2ν)], μ = E/[2(1 + ν)], where E is Young's modulus), εi are the principal strains, θ = ε1 + ε2 + ε3, and l, m and n are the Murnaghan third-order elastic constants. [17] Rotating equations (2) into principal strain coordinates, then and ε3 = 0 (plane strain). (The ± is + for the first subscript, − for the second subscript.) In this coordinate system, thermoelastic and poroelastic relative wave speed changes can be expressed as and where Vij0 are the wave speeds in the unstressed body, and the expressions are accurate for small ΔV/V0. Using the elastic strains of equation (9) results in identical expressions except with Ae(t) replacing A(t) and ky replacing (1 − ky). [18] Equation (14a) expresses the P wave component of the relative change in elastic wave speeds, whereas equations (14b) and (14c) express the S-wave components. Since the 1–2 axes do not (in general) line up with the x − y axes, only quasi-P, quasi-SV and quasi-SH waves exist in the x − y − z coordinate system. For the same reason, only quasi-Rayleigh and quasi-Love waves exist, and there remains coupling between them [e.g., Anderson, 1961]. However, approximate expressions for these anisotropic Rayleigh and Love wave components can be computed [Backus, 1970; Smith and Dahlen, 1973; Crampin, 1981]. It may be of use to note, also, that the anisotropic expressions of equations (14a), (14b), and (14c) share an isotropic component (terms multiplying sin kx) in addition to the anisotropic component (terms multiplying 1 − ky), so that there exists an isotropic part to the thermoelastic wave speed variations. Finally, one should note that all expressions make use of "average" Murnaghan constants that are assumed to be uniform throughout the entire half-space, an assumption that may not be very realistic. 3. Comparison of Seasonal Observations With Modeled Annual Variations [19] In the following subsections, I use the theory presented in section 2 to estimate whether it is possible for the three models discussed (thermoelastic, poroelastic hydrologic, and direct elastic hydrologic response) to explain the displacements and wave speed variations observed in southern California. The general approach taken is to use estimates of model parameters that are realistic but generous in terms of producing the desired effect, and to compare these "maximum realistic" model results with single observations that are representative of the southern California region. For example, as model parameters I choose temperature and water table amplitudes that are among the highest observed in southern California. The reasoning behind such a choice is that if these maximum realistic models are still unable to produce the desired effect, these models can be falsified more convincingly as being large contributors to the observed effect. It should be noted that the ability of some of these maximum realistic models to reach observed levels only implies that those models are potentially credible, not that they necessarily explain the observations everywhere throughout southern California (or elsewhere). Finally, as described in more detail below, model uncertainties are significantly higher for seismic wave speed variations than for displacements; nonetheless, as will be shown, all of the models require parameter choices outside the likely realistic range in order to produce wave speed variations that are in the observational range. GPS Observations [20] GPS displacement time series are well known to have seasonal variability [Dong et al., 2002] and a subsequent study by Prawirodirdjo et al. [2006] suggested that thermoelastic strains are a significant contributor to this variability. However, despite successfully matching much of the observed signal, Prawirodirdjo et al. [2006] do not attempt to quantitatively compare the observations with the thermoelastic model. In particular, they compare horizontal GPS displacements with calculated horizontal strains by arbitrarily scaling the amplitudes, without attempting to account for the fact that local displacements could have been determined from the modeled strain field. Comparison of equation (2a) with equation (11b) shows that ux has a spatial variability exactly out of phase with εxx. Furthermore, once k and Δt are determined, the amplitudes of the thermoelastic displacements are no longer arbitrary and can be calculated using equations (11a) and (11b). To make a quantitative comparison of the models with the observations, I apply the results of section 2.4 to the observations of Prawirodirdjo et al. [2006]. It should be noted that the direct elastic displacements from water table variability have been modeled previously by Bettinelli et al. [2008], and that Bawden et al. [2001] and Watson et al. [2002] have previously suggested groundwater seasonality to be important. 3.1.1. Thermoelastic Model Comparison [21] In using equation (11b) to calculate the annual variability in horizontal displacements, I choose representative values of parameters as before. To summarize, I take ν = 0.3, αth = 10−5°C−1, κ = 10−6 m2/s, ω = 2 · 10−7 s−1. To approximately match the data of Prawirodirdjo et al. [2006], I further take T0 = 10°C, k = 2π/(20 km), yb = 0.5 m. The value of T0 is chosen to match the high end of observed temperature variations in southern California (specifically, matching the annual Fourier component of the temperature record from Palm Springs as shown in Figure 3a); the value of k is chosen to fit the spatial variability in the observations; and the value of yb is chosen to fit the phase delay of the thermoelastic variations relative to the observed temperature variations, as described in equation (4). Here, and below, best-estimate values are provided. Figure 3Open in figure viewerPowerPoint Comparison of observed (solid blue) and modeled (dashed red) (a) T, (b) ux, and (c) ΔV/V0 as a function of time. All model results are shown at the horizontal position, x, for which the signal is largest. (Figure 3a) Smoothed temperature record from Palm Springs Airport (PSP) and modeled temperature with T0 = 10°C. (Figure 3b) Approximate averaged GPS N-S ground displacement from the 29 Palms region of Prawirodirdjo et al. [2006] and the modeled displacement of equation (15). The chosen GPS record is representative of displacements observed in the region; details regarding the GPS observation and site conditions are given by Prawirodirdjo et al. [2006]. (Figure 3c) Approximate averaged ΔV/V0 observed by Meier et al. [2010] in the Los Angeles basin and the modeled ΔVSV/VSV0 of equation (18) with m/μ = 104. [22] The horizontal displacement, ux, is evaluated at y ≈ 0 (close to the surface, where the GPS station is assumed to be anchored). Substituting these numbers into equation (11b) results in The observed amplitudes of GPS annual variability in the study region are about 2 mm. Thus, the thermoelastic displacements calculated could potentially represent a significant fraction (≈25%) of the observed displacements (see Figure 3b). Some of the parameter values (e.g., thermal diffusivity or thermal expansion coefficient) have some uncertainty and heterogeneity, and a slightly higher fraction of the observed signal could potentially be explained with appropriate modification of those values. One may note that after substitution of equation (5), equation (11b) is independent of k, has an amplitude that decays exponentially with yb, and is proportional to T0. It is likely that some fraction of the temperature fluctuation is recorded uniformly by the subsurface (resulting in a uniform temperature field without thermoelastic strains), implying a lower value of T0 and therefore explaining a smaller fraction of the observed displacements. Moreover, as stated above, the value of T0 used is already on the high end of observed values. It is therefore likely that the thermoelastic model achieves less than 25% of the observed displacement amplitude. [23] The spatial dependence of ux on −cos kx given a temperature variation with dependence on sin kx suggests that if thermoelastic variations are a significant component of the observed displacements then the maximum subsurface temperature variations occur spatially out of phase with the maximum observed displacements. This conclusion is counter to that assumed by Prawirodirdjo et al. [2006]. 3.1.2. Hydrologic Model Comparisons [24] For variations in water table level, I use equation (11b) to calculate the poroelastic horizontal displacement and the equivalently integrated form of equation (9a) to calculate the direct elastic horizontal displacement. As before, for parameter values I use E = 1.6 · 1010 Pa, α = 0.7 and ϕ = 0.15 to approximate sandstone values from Detournay and Cheng [1993]; and I use p0 = 2.9 · 104 Pa, equivalent to height of water of h0 = 3 m, to match observed water table variation in the (Los Angeles) region (see Figure 4a). Values of κhy for sandstone vary substantially, from 0.005 m2/s to 1.6 m2/s [Detournay and Cheng, 1993]; here, I use 0.4 m2/s as a nominal value, but note that this could be significantly different from the true value. I also note that the value of p0 (or equivalently, h0) also varies spatially in a significant manner; for example, some parts of southern California, such as the Mojave region, receive less rain than the Los Angeles region and have smaller values of p0 and h0. Figure 4Open in figure viewerPowerPoint Comparison of observed (solid blue) and modeled (dashed red) quantities from water table variability. (a) Water table record from California state well 04S12W36J001S in Rossmoor, California (http://www.water.ca.gov) and modeled water table with h0 = 3 m, corresponding to p0 = 2.9 · 104 Pa. (b and c) Comparisons equivalent to those in Figures 3b and 3c but for poroelastic effects. (d and e) Comparisons equivalent to those in Figures 3b and 3c but for direct elastic effects. Both Figures 4c and 4e use the less extreme choice of m/μ = 2 · 103. [25] Evaluating the modeled ux as previously, I obtain the results as plotted in Figures 4b and 4d for the poroelastic and direct elastic effects, respectively. As shown, both of these models have order-of-magnitude agreement with the amplitude of the observations, with a slight overprediction (by about 10%–30%). The phase of both models is reasonable but does not match the observations perfectly, partly due to the fact that the observations have subyearly Fourier components whereas the models have a single prescribed annual component. It should be noted that some of the overprediction is likely due to using a value of p0 larger than is appropriate for the 29 Palms region where the GPS observations are taken; some of the overprediction may also be due to the fact that the direct elastic and poroelastic effects tend to cancel each other (they would have exactly opposite signs for large values of κhy), with the true effect being equivalent to the difference in the two signals. Finally, the overprediction may also be because the 29 Palms observations are not entirely representative of displacements in the southern California region; for example, for a similar region, Watson et al. [2002] find peak-to-peak displacements on the order of 6 mm, rather than the 4 mm of Prawirodirdjo et al. [2006]. 3.1.3. Summary of GPS Displacement Comparison [26] Based on the above considerations, I conclude that thermoelastic variations may be responsible for an observable fraction of the annual variability of horizontal GPS displacements but that it is not likely to explain the entire annual signal even in places with large temperature fluctuations. I further conclude that in regions with significant water table variability, those variations could potentially explain the full signal through either the