A viscoelastic model for groundwater level changes in the Cho-Shui
River alluvial fan after the Chi-Chi earthquake in Taiwan
Yun-Bin Lin,1,2 Yih-Chi Tan,3Tian-Chyi Jim Yeh,4Chen-Wuing Liu,3 and Chu-Hui Chen5
Received 20 June 2003; revised 21 February 2004; accepted 2 March 2004; published 29 April 2004.
[1] A viscoelastic model is developed to simulate the groundwater level changes in the Cho-Shui River alluvial fan in Taiwan after the Chi-Chi earthquake. An analytical solution is derived with the assumption that no leakage occurred in confined aquifers during the coseismic period. The solution is used to analyze the data collected from a high-density network of hydrologic monitoring wells in the Cho-Shui River alluvial fan. The simulated groundwater level changes agree with the observations. The viscosity coefficient of the model was found to correlate with the hydraulic conductivity of the aquifer. The field observations and the simulations reveal the influence of geological structures and heterogeneity on the groundwater changes and locations of sediment liquefactions in the alluvial fan during the Chi-Chi earthquake. Possible applications to imaging subsurface hydraulic heterogeneity are discussed using information about groundwater level changes induced by earthquakes. INDEXTERMS: 1829 Hydrology: Groundwater hydrology; 1878 Hydrology: Water/energy interactions; 7260 Seismology: Theory and modeling; KEYWORDS: Chi-Chi earthquake, groundwater level changes in wells, viscoelastic model, delineation, sand/gravel formation
Citation: Lin, Y.-B., Y.-C. Tan, T.-C. J. Yeh, C.-W. Liu, and C.-H. Chen (2004), A viscoelastic model for groundwater level changes in the Cho-Shui River alluvial fan after the Chi-Chi earthquake in Taiwan, Water Resour. Res., 40, W04213,
doi:10.1029/2003WR002412.
1. Introduction
[2] At 01:47012.600on 21 September 1999, an earthquake
of magnitude 7.3 on the Richter scale (ML) occurred in the
middle of Taiwan, caused by the motion of the Chelungpu thrust fault. The epicenter was located at 23.87N and 120.75E, near the town of Chi-Chi, which is about 10 km away from the Cho-Shui River alluvial fan [Ma et al., 1999; Yu et al., 2001]. A second earthquake of a magnitude of 6.9 took place at 19:52050.000 on 26 September 1999 with the epicenter at 23.86N and 121.00E. The Cho-Shui River alluvial fan is located in the midwestern plan of Taiwan and has an area about 1700 km2. A series of west vergent thrust faults and anticlines have been mapped near the eastern side of the fan (Figure 1). The Cho-Shui River flows from east to west through two anticlines, the Dulliu and Baqua hills [Lin et al., 2000, 2001], and is the main recharge area of the fan. A groundwater monitoring network system (GMNS) was established and has been maintained across fan since 1992. The network consists of 70 evenly distributed hydrological stations, where a total of 188 monitoring wells are installed at various depths from 24 to 306 m below ground surface.
Each well is screened at only one depth, and the water levels at all of the wells have been recorded automatically at every hour [Hsu, 1998]. The records of the GMNS showed the groundwater level rises ranged from 1 m to 7 m in some wells and falls from2 m to 11 m in the others after the Chi-Chi earthquake [Water Conservancy Agency, 1999a].
[3] During the past few decades, efforts have been made
to identify the relation between field strain and groundwater level changes caused by earthquakes. Numerous studies [Roeloffs, 1996] have suggested that the poroelastic theorem by Biot [1956a, 1956b] can be used to describe the local pore pressure change induced by earthquakes. Many attempts also have been made to use the groundwater level change as a precursor of the earthquake [Ohno and Wakita, 1997]. However, because of insufficient monitoring wells [Roeloffs, 1998] and the scale disparity between models and the aquifer [Grecksch et al., 1999], the relation between earthquakes and changes in groundwater levels remains uncertain.
[4] The well-instrumented well fields in the fan, and the
vast amount of data recorded by the GMNS during the Chi-Chi earthquake provide a unique opportunity to study the relation between earthquakes and groundwater level changes. The objective of this study is to construct a mathematical model to describe and understand the ground-water level changes after the Chi-Chi earthquake in the Cho-Shui River alluvial fan. Specifically, we develop a viscoelastic model based on the local hydrogeological condition of the alluvial fan, derive an analytical solution of the model, and use it with field data to simulate the evolution of the groundwater level changes after the earth-quakes. Finally, we discuss the importance of geological structures, and heterogeneity on the observed groundwater
1
Council of Agriculture, Executive Yuan, Taipei, Taiwan.
2Also at Department of Bioenvironmental Systems Engineering,
National Taiwan University, Taipei, Taiwan.
3Department of Bioenvironmental Systems Engineering and Hydrotech
Research Institute, National Taiwan University, Taipei, Taiwan.
4Department of Hydrology and Water Resources, University of Arizona,
Tucson, Arizona, USA.
5Department of Civil Engineering, ChungKuo Institute of Technology,
Taipei, Taiwan.
Copyright 2004 by the American Geophysical Union. 0043-1397/04/2003WR002412$09.00
level changes and sediment liquefaction during the earth-quakes. Possible applications of the model to identification of aquifer heterogeneity are also explored.
2. Theory of Propagation of Elastic Waves in Aquifers
[5] Hassanizadeh and Gray [1979a, 1979b, 1980]
adop-ted the method proposed by Coleman and Noll [1963] to linearize the momentum equations of two-phase flow by using the entropy inequality of thermodynamics. Linearized momentum equations of the liquid phase and the solid phase can be expressed as 1 ef rs@ 2us k @t2 1 ef rsgs kþ 1 ef rkpf sskl;lþ R fud k¼ 0; ð1aÞ and efrf @2uf k @t2 efr fgf kþ efrkpf Rfudk ¼ 0 ð1bÞ
respectively, with the assumptions that (1) there is no phase change; (2) the fluid phase is macroscopically nonviscid; (3) the temperature gradient is considered only as heat conduction; (4) the acceleration of the convection term is neglected; and (5) the capillary effect is negligible. In equations (1a) and (1b),eais the void fraction with respect
to thea phase, f denotes the liquid phase, s denotes the solid phase,ra represents intrinsic mass density of thea phase
[M/L3], uak is the displacement vector of thea phase in the
k direction [L], g is the gravity [L/T2],rkp f
denotes pressure gradient of the liquid phase in the k direction [M/(L3T2)],skl,ls
denotes the stress tensor gradient with respect to the solid phase in l direction [M/L3T2)], Rf is a linear coefficient tensor describing the momentum exchange due to the relative velocity of the liquid phase and the solid phase [M/(L3T)]. The coefficient tensor becomes a scalar when-ever the reference coordinates are parallel to the principle axes of the equations. The term ukdis the relative velocity
vector of the liquid and solid phases in the k direction [L/T]. [6] When equations (1a) and (1b) are used to describe the
propagation of elastic waves in fluid-saturated porous media, Biot [1956a, 1956b], Berryman et al. [1988], and Jeng and Lee [2001] showed that an extra term is needed to fully represent the dynamic behavior of the media. The term is related to the relative acceleration of the liquid and solid phases. After adding this term, the momentum equations become 1 ef rs@ 2us k @t2 1 ef rsgs kþ 1 ef rkpf sskl;l þ Rfud kþ A sad k ¼ 0; ð2aÞ and efrf @2uf k @t2 efr fgf kþ efrkpf Rfudk A sad k¼ 0: ð2bÞ
In which As is an additional apparent mass density of the solid phase [M/L3]; and akdrepresents the relative
accelera-Figure 1. The geography, geology, including faults and epicenters, and locations of the well stations in the study area.
tion of the liquid phase to the solid phase in the k direction [L/T2]. Equations (2a) and (2b) can be compared with the coupled momentum equations for solids and fluids of Biot’s theory of propagation of elastic waves in a fluid-saturated porous medium: @2 @t2 r11u s kþ r12u f k þ 1 ef rkpf sskl;l bu d k ¼ 0; ð3aÞ and @2 @t2 r12u s kþ r22u f k þ efrkpfþ budk¼ 0 ð3bÞ where the coefficient b is related to the intrinsic perme-ability, k[L2], the fluid viscosity, m [M/(LT)], and the porosity as b =me2f/k. Here the gravity term is also omitted
if a horizontal plane model is considered. The comparison yields that (1)r12= As,r11=rses Asandr22=rfef As;
(2) b = Rf; (3) r11 + r12 is thus equal to the bulk density
of the solid phase per unit volume of porous media, esrs;
(4)r12+r22, is equivalent to the bulk density of the liquid
phase per unit volume of porous media, ef rf; and (5) the
term, akd, is equal to @
2
@t2 (ukf uks).
[7] If equations (3a) and (3b) are isotropic functions, the
representation theorem [Wang, 1970a, 1970b; Drew and Passman, 1999] explains the necessity of the term, As. On the basis of the representation theorem, an isotropic function should include all the objective variables. Variables can be regarded as objective only when their values are frame-indifferent. Because akdis an objective variable, Asakdshould
be included in (2a) and (2b) such that the isotropic property of the equations can be fully specified. Physically, the term,r12, can be regarded as an additional mass, induced
by the oscillation of solid particles in the liquid phase. A similar theorem explaining an apparent mass induced by a rigid body oscillating in the fluid can also be found in the potential flow theorem of the fluid mechanics [Lamb, 1945]. [8] Using the Lagrangian formulation, Berryman et al.
[1988] developed a complete theory of elastic wave propa-gation in partially saturated porous media. Equations of motion were derived for solid, liquid, and gaseous constit-uents, including the effects of their interactions. Assuming the effects of capillary pressure changes are negligible, Berryman et al. [1988] showed that the equations could be simplified to a form similar to Biot’s equations. Using the Helmholtz transformation, they decoupled the equations into
r2c¼k 2 s w2 @2c @t2 ð4aÞ
for one distortional wave and
r2 A¼ k2 w2 @2A @t2 ð4bÞ
for two dilatational waves. In equations (4a) and (4b), A+is
a scalar potential for the solid phase with the faster wave speed, and Ais a scalar potential for the liquid phase with the slower wave speed; ks, and k±are wave vectors [1/L];c
is a vector potential for the solid phase; and w is the vibration frequency [1/T].
[9] If we assume equations (3a) and (3b) are isotropic and
take the divergence of the equations, we have
r2 ss 1 ef r2 pfþ rpf re f ¼ @2 @t2ðr11eþ r12VÞ b@ @tðV eÞ u dr b; ð5aÞ and efr2pf rpf ref ¼ @2 @t2ðr12eþ r22VÞ þ b @ @tðV eÞ þ udr b; ð5bÞ where e = r uk s andV = r uk f
. According to the Biot theorem, the pore pressure obeys the diffusion equation in the poroelastic model. Equations (5a) and (5b) describe the fully dynamic behaviors of the solid and the liquid phases. The summation of equations (5a) and (5b) leads to
r2 ss e spf efpf ¼ @ 2 @t2½ðr11þ r12Þe þ rð 22þ r12ÞV : ð6Þ Using the mixture theorem [Green and Naghdi, 1965], we can express equation (6) as two wave equations:
r2 ss espf ¼ @ 2 @t2½ðr11þ r12Þe ¼ r s@2 @t2ðeseÞ ffi rs@2 @t2 r esu s ð Þ ½ ð7aÞ r2 e fpf ¼ @ 2 @t2½ðr22þ r12ÞV ¼ rf @2 @t2 efV ffi rf @ 2 @t2 r efuf : ð7bÞ
This theorem implicitly assumes that terms such as rpf ref,@
2
@t2(r12e +r12V) and udr b are small compared with
the other terms in equations (5a) and (5b), and they can be neglected in equations (5a) and (5b) during the coseismic motion period. Because of the omission of the term,
@2
@t2(r12e +r12V), equations (7a) and (7b) are quasi-dynamic
descriptions [Jeng and Lee, 2001] of equations (5a) and (5b). Additionally, the b@t@(V e) term is omitted in the summation but its effect can be included using a constitutive equation to be discussed later.
[10] Prevost [1982] suggested that a single-phase
descrip-tion of porous media behavior is adequate when the loading rate is much higher than the diffusion rate. This implies that if no efflux of the pore fluid from porous media takes place during the coseismic motion period and the fluid is com-pressible, then only equation (7b) is needed for describing the quasi-dynamics of wave propagation in porous media. As a result, we assume that no leakage occurs from one confined aquifer to another in the fan during the coseismic compression. This is supported by hydrogeological data supplied by the Central Geological Survey [1994, 1999], shown in Figure 2, and we can use only equation (7b) in our following analysis.
[11] Suppose that the liquid phase in our study area obeys
the Maxwell-Fluid model (see Figure 2) [Flu¨gge, 1975], a constitutive relationship between the stress and the strain
can be used for describing the pore pressure dissipation caused by the relative velocity of the liquid phase and the solid phase: efpf þh C @ @t efp f ¼ h@ @t r efu f : ð8Þ
In equation (8), the viscosity coefficient, h [M/(LT)], is defined as the ratio of the pressure to the strain rate and C is
the bulk modulus of the liquid phase [M/(LT2)]. The
advantage of using equation (8) is that knowledge of the exact motion of the solid phase is not necessary, which is implicitly included in the viscosity coefficient in equation (8). The coefficient can be related to well-defined parameters, such as hydraulic conductivity.
[12] For the one-dimensional model, the horizontal
dis-placement of the liquid phase per unit volume of the porous medium can be expressed as
efuf ¼ xegxei lxwtð Þ; ð9Þ
where x is a constant [L]; g is an attenuation coefficient [1/L]; x is the distance [L]; l is the wave vector, [1/L], and w is the vibration frequency [1/T]. The strain of the liquid phase per unit volume of porous media then is
r efuf
¼ g þ lið Þ efuf
: ð10Þ
The strain rate of the liquid phase per unit volume of porous media therefore becomes
@r efuf @t ¼ iw efu f g þ li ð Þ: ð11Þ
In addition, we express the relation between the strain rate of the liquid phase per unit porous media and the stress of the liquid phase per unit porous media as
L@r efu f
@t ¼ efp
f: ð12Þ
Figure 2. Hydrogeological profile of the alluvial fan and the derived conceptual model.
Substituting this relation into equation (8), we obtain
L¼ h
1whi C
: ð13Þ
According to equations (12) and (13), the phase of the pore pressure per unit volume of porous media has a lag of q = tan1 wh C behind the phase of the strain rate of the liquid phase per unit volume of porous media. Upon substituting equations (10), (11), (12) and (13) into (7b), equation (7b) becomes whi 1whi C r2 efuf ¼ rf @2 @t2 efu f : ð14Þ
Comparing equation (14) with equation (4b), we can set A=efufto get k2= w
2
v2 f
(1 + Qi) where Q =whC. From equation (9), we can show that
@2 @t2 efu f ¼ iw@ @t efu f : ð15Þ
Using equation (15) in equation (14), equation (14) can be expressed as @2A @x2 rf h @A @t 1 v2 f @2A @t2 ¼ 0; ð16Þ where vf = ffiffiffiffiffiffiffiffiffiffi C=rf p
is the velocity of the wave propagation in the liquid phase. Depending on the equations (11), (12), and (13), the pressure per unit volume of the fluid can be expressed in a similar form as equation (16). This is our viscoelastic model to be used in the following analysis. The viscoelastic model considers the inertial effect that is neglected in the poroelastic model. Eringen [1980] and Bardet [1992] showed that a viscoelastic model, instead of the poroelastic model, could describe the dynamic response of the saturated porous media.
3. Hydrogeology of the Study Area
[13] The Cho-Shui River’s alluvial fan consists of several
layers of Holocene to Pleistocene sands and gravels that formed the three aquifers separated by marine mud [Central Geological Survey, 1994, 1999]. The rising and falling mean sea levels caused by global climate change late in the Quaternary Period formed the layered structure of the alluvial fan. Massive gravels, which comprise many layers of the upper fan, tend to thin out toward the west of the fan, while the mud layers thicken.
[14] The Central Geological Survey [1994, 1999]
con-structed 12 hydrogeological profiles for the alluvial fan; three aquifers (aquifers 1, 2, and 3) were identified in the analysis of the profiles (Figure 2a). Aquifer 2 can be divided horizontally into a confined part on the tail fan, a partially confined part on the mid fan and an unconfined part on the upper fan. Gravelly sediments were found on the middle-upper fan and in the partially confined and unconfined parts of the aquifers, while arenaceous sediments were found only
in the confined parts of the aquifer. Furthermore, an interface was identified between the gravels and the arenaceous sediments. This interface separates the partially confined part from the confined part of the aquifer. The confined portion pinches out below the shoreline (Figure 2b).
[15] Rupturing of the Chelungpu fault in the Chi-Chi
earthquake led to several meters of oblique thrust of the hanging wall relative to the footwall. According to GPS measurements, only minor subsidence (<0.5 m) was ob-served in the Cho-Shui River alluvial fan after the Chi-Chi earthquake [Department of Land Administration, 2000]. The groundwater level in the Cho-Shui river alluvial fan, however, changed as much as ±10 m right after the Chi-Chi earthquake [Chi-Chia et al., 2000; Water Conservancy Agency, 2000]. The spatial distribution of general ground-water level changes (i.e., rises and declines) between 01:00 and 02:00 on Sep. 21 are shown in Figure 3. On the basis of Figure 3, negative water level changes were recorded in wells located in a narrow, belted area between the Chelungpu and Changhwa thrust faults. This volumetric expansion zone in the footwall perhaps is caused by the dragging effect of thrusting. Wells that recorded positive water level changes are generally located further away from the ruptured fault. The magnitude of the rise of the water level tends to increase from the upper fan toward the mid fan, and then decrease toward the tail fan [Chia et al., 2001; Wang et al., 2001; Lee et al., 2002]. In wells adjacent to the western side of the Changwha thrust fault, the groundwater levels declined in one aquifer and rose in all the other aquifers. Such a distribution of groundwater level changes was found inconsistent with the calculated strain fields reported by Ma et al. [2001] and Huang [2000]. The inconsistency appears to suggest that ground-water levels may also be affected by the hydrological heterogeneity in the aquifers in the fan in addition to the forces of the earthquake.
4. Conceptual Model for the Study Area
[16] The Chelungpu fault ruptured from the south to the
north during the Chi-Chi earthquake and induced pressure waves, propagating from east to the west and almost parallel to the contour lines of the groundwater levels recorded before the Chi-Chi earthquake. Therefore the 85 km long rupture zone that formed during the Chi-Chi earthquake supplies a western-eastern direction pressure pulse simulta-neously to the fan (Figure 3). Given the timescale of the model (a daily average) and the dozens of second coseismic periods, the pressure induced by the earthquake is treated as an impulse pressure. Suppose that no leakage took place from the confined part to other aquifers during the coseis-mic compression. Since the thickness of the aquifer decreases from east to west, conservation of momentum suggests that when the wave propagates from east to west, the change of the velocity will be significant toward the pinch out of the aquifer. Thus the inertial influence on the groundwater level fluctuations increases. The inertial term must be considered at the western side of the interface; our viscoelastic model, equation (16), thus will apply. Con-versely, on the basis of the physical properties of gravelly aquifers and the hydrogeological structure of the aquifers, the fluid behavior after the earthquake at the eastern side of the interface is assumed to obey the diffusion equation.
[17] To use equation (16) to describe the motion of the
fluid phase within the arenaceous confined aquifer, the boundary conditions of equation (16) are set to be a zero strain rate at the pinch out (x = L, the distance between the interface and the pinch out) and an impulse pressure at the interface (x = 0) (see Figure 2c). The initial conditions for equation (16) are
Aðx; 0Þ ¼ @
@tAðx; 0Þ ¼ 0: ð17Þ Equation (17) implies that the velocity and displacement of the liquid phase are zero before the earthquake.
[18] Taking the Laplace transformation of equation (16)
to project t into s and A(x, t) into A(x, s), and then using
the initial conditions, equation (17) becomes
A 00 x; s ð Þ r f hsAðx; sÞ 1 v2 f s2A ðx; sÞ ¼ 0: ð18Þ
We assume the solution takes the form
A¼ M sð Þelx: ð19Þ
After substituting equation (19) into equation (18), we have
l¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rf hs 1þ h Cs s ¼ 1 vf ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s C hþ s s : ð20Þ
Because when x! 1, A! 0, the negative l value is the
correct solution.
[19] When the impulse pressure Pd(t) is applied at x = 0
and t = 0, the strain at x = 0 is
A0ð0; tÞ ¼ Pd tð Þ h Cþ t
h : ð21Þ
Figure 3. The distribution of the rising (solid circles), the falling (squares), and the partially rising (open circles) groundwater levels of well stations and the groundwater level before the Chi-Chi earthquake in aquifer 1.
The Laplace transform of equation (21) yields A 0 0; s ð Þ ¼ P1 C¼ lM sð Þ: ð22Þ
After M(s) is determined using the boundary condition, it is substituted into equation (19). We then have
Aðx; sÞ ¼ Pvf C ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s C hþ s s e x vf ffiffiffiffiffiffiffiffiffiffiffi sC hþs ð Þ p : ð23Þ
The Laplace transformation of f(t) = I0(v
ffiffiffiffiffiffiffiffiffiffiffiffiffiffi t2 x2 p )H(t x) is f (s) = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1 sv ð Þ sþvð Þ p ex ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðsvÞ sþvð Þ p
[Farrell and Ross, 1971]. Meanwhile, the translation property projects f (t) = I0(v ffiffiffiffiffiffiffiffiffiffiffiffiffiffi t2 x2 p )evtH(t x) into f (s) = ffiffiffiffiffiffiffiffiffiffiffiffiffi1 s sþ2vð Þ p ex ffiffiffiffiffiffiffiffiffiffiffiffiffis sþ2vð Þ p . Usingv = C 2h and x = x
vf in equation (18), the transverse
Laplace transformation of equation (23) can be expressed as
Aðx; tÞ ¼ Pvf C I0 C 2h ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t2 x vf 2 s 0 @ 1 AeC 2htH t x vf : ð24Þ
The pore pressure per unit volume of porous media can be obtained by applying equation (24) to equation (11) and expressed as efpf ¼ Pebt 1þ 1 Q 2! x vf tI1ð ÞQ 2D3=2 CI1ð ÞQ 4hD1=2 C 8h t I½0ð Þ þ IQ 2ð ÞQ D H t x vf þ Z; ð25Þ
where D = t2 (x/vf)2; H( ) is the Heaviside function;Q = b
ffiffiffiffi D p
; b = C/(2h), and In( ) is the nth-order Bessel function.
[20] The constant Z introduced by the phase lag in
equation (25) can be obtained by applying the boundary condition of zero strain rate, that is, @
@t(r A(L, t)) =
efpf
L = 0. Therefore the pore pressure per unit volume
of porous media induced by unit pressure, Unitf(x, t) when t > (x/vf), is given as Unitfðx; tÞ ¼ e bt 1þ 1 Q 2! tI1ð ÞQ 2D3=2 CI1ð ÞQ 4hD1=2 C 8h t I½0ð Þ þ IQ 2ð ÞQ D ( ) x vf tI1ðQLÞ 2D3=2L CI1ðQLÞ 4hD1=2L C 8h t I½0ðQLÞ þ I2ðQLÞ DL ( ) L vf 8 > > > > > < > > > > > : 9 > > > > > = > > > > > ; H t x vf ð26Þ where DL= t2 (L/vf)2andQL=b ffiffiffiffiffiffi DL p . From equation (26), whenh is smaller, the pore pressure dissipates faster.
[21] In order to apply the solution, equation (26), to the
study area, the pressure impulse at the interface between the
gravels and sands must be specified. To resolve this issue, we use the following approach. The force F of the earth-quake leads to a sudden increase in pore pressure at the interface. On the basis of Biot’s derivation [Biot, 1941; Ge and Stover, 2000] and our knowledge of the geohydrology of the fan, whenever the force F is removed, the dissipation of the pore pressure, efpbc
f
, on the eastern side of the interface will obey the diffusion equation,
@ efp f bc @t ¼ K Ss @2 e fp f bc @x2 þ F Bd xð Þd t t 0 ð Þ; ð27Þ
where Ss is the specific storage [1/L]; K is the hydraulic
conductivity of the aquifer [L/T], and B is the thickness of the aquifer. If the initial condition is given as pbcf(x, 0) = 0,
the corresponding solution for an unbounded medium is
pfbcðx; tÞ ¼ F Bexp x2S s 4Kt 2ef ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pt K=Sð sÞ p : ð28Þ
When the impulse pressure is applied to the liquid phase, the displacement takes place simultaneously along the direction of the force, and the energy associated with the impulse is transmitted via the wave propagation from the upper fan to the tail fan. The geological structures and heterogeneity caused the high pore pressure at the interface in the aquifer. Because of the unbalanced pressure head, the liquid phase flows from high to low pressure. After the impulse pressure was removed, i.e., after the Chi-Chi earthquake, groundwater flowed from the sand and gravel interface in an upstream and downstream direction in the Cho-Shui River alluvial fan. The energy is proportional to square of the displacement; the displacement will decrease with the time and distance of the wave propagation because of the viscosity between the liquid phase and the solid grains. If the upstream flow obeys the diffusion equation, then the change in the pore pressure at the interface, x = 0, is pfbcð0; tÞ ¼ 1 2ef ffiffiffiffiffi pt p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiF=B K=Ss ð Þ p : ð29Þ
The pore pressure in the arenaceous sediments, pres
f
, therefore can be expressed as the convolution form of pbc f and Unitf as pfresðx; tÞ ¼ Z t 0 pfbcð0; tÞUnitf x; t t ð Þdt: ð30Þ
5. Results and Discussion
[22] To apply the model to the study area, observed
groundwater levels from four sets of wells are used (Figure 4). Each well set consists of one well on the east side of the interface and another on the west side. Each well is labeled using a number followed by a letter and another number (e.g., 3W1). The first number denotes the aquifer where the well is screened (i.e., 3 for aquifer 3); the second
letter indicates the location of the well in relation to the interface (i.e., W for west of the interface), and the last number is the set number of the well (i.e., 1 for the first set). [23] To simulate the dissipation of the pore pressure in the
arenaceous sediments, first the pore pressure dissipation at the interface (i.e., F/b) is found by fitting equation (28) to the groundwater level records of the eastern well with known Le, K and Ss. The variable, Le, denotes the horizontal
distance in the western-eastern direction between the inter-face and the well. The values of K and Ssare obtained from
aquifer tests at the wells. Once the pore pressure dissipation at the interface is known and the east-west horizontal distance between the interface and the well, Lw, is given,
the h value in equation (16) is determined by fitting
equation (30) to the observed groundwater level of the western side well. The bulk modulus of the fluid phase,
C, is set to be 2.25 109Pascal, and the distance from the
interface to the pinch out, L, is set to be 30,000 m. [24] According to the daily average records, the largest
groundwater level changes in the wells were located at the interface between the gravels and the sands. To test the model, two well sets were chosen, which include 3E1, 3W1, 3E2, and 3W2, shown in Figure 4c. The records of wells 3E1 and 3E2 were fitted by equation (28), and the records of wells 3W1 and 3W2 were fitted using equation (30). Table 1 lists the best fit parameters’ values and associated input parameter values. The simulated and observed groundwater levels as a function of time are plotted in Figures 5 and 6. In general, the agreements are excellent; the small differences between the observed records and the simulated results may be attributed to the simplified one-dimensional model, as well as the no-leakage assumption. Figure 4. Changes in groundwater level, hydraulic conductivity distribution, and the locations of the
interface, pinch out, and observed liquefaction areas.
Records of another two sets of wells on the southern part of the fan were also analyzed, shown in Figures 4b and 4c as 3E3, 3W3, 2E1 and 2W1. Wells 3W3 and 2W1 are located at the same horizontal position, but with different screen depths. Simulated results for wells 3E3, 3W3, 2E1 and 2W1 are plotted against observed data (Figure 7). The observed dissipation of the pore pressure at well 2E1 is faster than simulated. This discrepancy may be attributed to the mul-tidimensional diffusion process and leakage between aqui-fers, which were neglected in the one-dimensional model. Well 2E1 is located near location X, which is a groundwater discharge area of the Cho-Shui River alluvial fan (Figure 4c). A sudden increase of flow in the channel of the Cho-Shui River, downstream from X, was observed after the Chi-Chi earthquake [Water Conservancy Agency, 1999b]. According to the Water Conservancy Agency [2000] report and Lee et al. [2002], the increase in river discharge was caused by groundwater discharge from aqui-fer 1, resulting from vertical water movement from aquiaqui-fer 2 to aquifer 1. Since no aquitard exists between aquifer 2 and aquifer 1, the vertical water movement is significant. As a result, the reduction in the observed head of aquifer 2 is much more rapid than the simulated water level at Well 2E1. The vertical water movement also may have occurred from aquifer 3, but its magnitude may have been reduced by the aquitard found between aquifers 2 and 3. The absence of the aquitard at Well 2E1 and its presence at Well 3E3 may provide a plausible explanation for the faster decrease of the
groundwater level observed at Well 2E1 (Figure 7), com-pared with the simulated result. There is no significant difference between the observed and simulated hydrographs at Well 3E1 (Figure 5). Other evidence of vertical water movement can be found at location Y in Figure 4c. At this location, an aquitard exists between aquifers 2 and 3. This aquitard’s eastern boundary is delineated in Figure 4c, which resembles the shape of the contour of the groundwa-ter level changes in aquifer. The similarity in shape may suggest up coning of the groundwater in the area. Therefore we conclude that the vertical water movement is significant in the unconfined aquifer of the fan.
[25] A plot of h/C values obtained from the calibration
and known Lw/K values is shown in Figure 8. These K
values are those estimated from aquifer tests at the wells on the western side of the interface. A best fit line is also shown, which is given as
h C ¼ 7:4 103 Lw K 0:6: ð31Þ
The correlation coefficient is 0.96. The viscosity dissipation is proportional to the distance and inversely proportional to the hydraulic conductivity. On the basis of this relation, we surmise that the further the well is located from the interface, the less the well is affected by the viscosity dissipation. If hydraulic conductivity increases, the viscosity
coefficient h decreases and the dissipation of the pore
pressure will occur quickly. Table 1. Fitting Parameters and Hydraulic Properties of the Wells Used in This Study
Well Leor Lw, m Positiona Specific Storage Ss, b 1/m Hydraulic Conductivity K ,bm/d Horizontal Distance Between the Eastern Well and the Western
Well, m hC, days X Component Y Component 2E1 0 185350 2624184 0.80e-4c 12.355 6,682 0.19 2W1 6600 178717 2624989 30.499 6,682 0.19 3E1 450 194052 2656100 1.50e-4 26.870 3,934 0.64 3W1 3500 190120 2656250 56.572 3,934 0.64 3E2 2000 196133 2649778 2.30e-4 52.445 8,613 0.00025 3W2 6000 187624 2648441 67.133 8,613 0.00025 3E3 680 191168 2627781 0.50e-4 29.894 12,746 7.2 3W3 12000 178717 2624989 11.568 12,746 7.2
aTransverse Mercator of 2, central meridian = 121E. b
Data were obtained by the well pumping tests [Water Conservancy Agency, 1999a].
cRead 0.80e-4 as 0.80
104.
Figure 5. The comparison of observed and simulated
groundwater levels in wells 3E1 and 3W1.
Figure 6. The comparison of observed and simulated
[26] Table 1 shows that the range of the fittedh/C value is
between 7.2 and 0.00025 days. This is indicative of visco-elastic behavior of the geological materials in the fan during the earthquakes. Assuming the bulk modulus of water, C, is 2.25 109 Pascal, the calculated viscosity coefficient, h, ranges between 4.86 1010Pa s and 1.40 1015Pa s. The
viscosity coefficient of water at 20C is 103 Pa s and the mantle of the earth has the viscosity coefficient of 1020Pa s. The large values of the calculated viscosity coefficient for the fan indicate that the material of the fan during the earthquakes could be in a stage of debris flow (perhaps, liquefaction). However, such a high viscous coefficient may be a result of our one-dimensional analysis with simplified assumptions.
[27] Notice that the Biot’s poroelastic model considers
fluid-saturated material with separate liquid and solid phases. Using the hydraulic conductivity and specific stor-age of the aquifer, the Biot’s poroelastic diffusion equation yielded less satisfactory results than the viscoelastic model with viscosity coefficient. Furthermore, because the geolog-ical structure (Figure 2) shows that the thickness of the aquifer becomes thinner from the upper fan to the pinch out, the inertial effect plays an important role on the momentum exchange of the fluid moving from the upper fan to the pinch out of the aquifer. Our viscoelastic model considers the inertial effect, which is neglected in the poroelastic model.
[28] After the Chi-Chi earthquake, liquefaction of
sedi-ments was observed at several locations in the alluvial fan (see Figure 4a). While many factors control the liquefaction, we postulate that this liquefaction was caused by the high pore pressures induced by the earthquakes at the interface, and the subsequent pressure propagation throughout the aquifers. The phenomenon is analogous to a choke in open
channel flow, or a water hammer in pipe flow. Several liquefactions of sediments were observed at the upper fan (i.e., A, B1, B2 and B3 in Figure 4a) where the aquifer is unconfined and strong vertical flow took place. Liquefac-tion of the sediment was also observed at the northern fan nearby the Taiwan Strait (see C and D in Figure 4a). At the northern fan, no impermeable clay sediments have been found on top of the aquifer along the coastline. The absence of clay sediments is possibly caused by wave erosion [Central Geological Survey, 1994, 1999]. The excess pore pressure caused by the earthquake thus was released at these outlets and induced the liquefaction. However, no sediment liquefaction was detected on the southern part of the fan, which is also nearby the coastline and much closer to the pinch out. This phenomenon appears to support our assumption of a zero strain rate at the pinch out. Most certainly, the interface-induced high pore pressure increases
Figure 7. The comparison of observed and simulated groundwater levels in wells 3E3, 3W3, 2E1, and
2W1.
Figure 8. The regression relation ofChandLw
K and the 95%
confidence interval of the regression.
the risk of sediment liquefaction, but the local hydrogeology ultimately controls the liquefaction potential.
[29] Finally, our analysis has demonstrated that the
prop-agation of groundwater waves induced by the earthquakes is controlled by the geological structures and hydraulic prop-erties of aquifers. The results of the analysis suggest that it is possible to use earthquake-induced groundwater changes to identify structures of groundwater basins and aquifer heterogeneity. Groundwater level changes caused by earth-quakes from different epicenters may provide data sets of naturally occurring hydraulic tomography surveys of the aquifers, surveys that are analogous to the hydraulic tomography proposed by Yeh and Liu [2000] and Liu et al. [2002].
6. Conclusions
[30] This study developed a viscoelastic model to
de-scribe the pore pressure dissipation after impulse pressures induced by two significant earthquakes. This model was used to reproduce the postearthquake groundwater level changes observed in monitoring wells in the Cho-Shui River alluvial fan. Once the direction of the wave propaga-tion was known, the viscosity coefficient of the model was estimated using the observed groundwater level changes in two monitoring wells in the direction of the wave propaga-tion, and the interface position of sand and gravel formation in the confined aquifer. Using this approach, the only data required is the hydraulic conductivity and the specific storage, both of which are obtained by field tests; knowl-edge of the magnitude of the earthquake is unnecessary. Results of the study show that the estimated viscosity coefficient is strongly correlated to the hydraulic conduc-tivity of the aquifer. Most certainly, local hydrogeological conditions influence liquefaction susceptibility. However, our model appears to suggest that the geological interface-induced high pore pressure may be the driving force. Finally, we conclude that geological heterogeneity plays a significant role in groundwater level changes and possibly sediment liquefactions induced by earthquakes. Our study may also catalyze the use of earthquake-induced ground-water level changes to delineate large-scale aquifer hetero-geneity and geological structures.
[31] Acknowledgments. This study was supported by the National Science Council of TAIWAN under contract NSC 89-2116-M-002-053. The authors are grateful to the Water Resources Agency and the Central Geological Survey of TAIWAN for providing the field data. The third author is funded by NSF and SERDP of USA by a grant EAR-0229717. The authors also extend their gratitude to Martha P. L. Whitaker for her editing of this document.
References
Bardet, J. P. (1992), A viscoelastic model for the dynamics behavior of saturated pore-elastic solids, J. Appl. Mech., 59, 128 – 135.
Berryman, J. G., L. Thigpen, and C. T. Chin (1988), Bulk elastic wave propagation in partially saturated porous solids, J. Acoust. Soc. Am., 4, 360 – 373.
Biot, M. A. (1941), General theory of three-dimensional consolidation, J. Appl. Mech., 12, 155 – 164.
Biot, M. A. (1956a), Theory of propagation of elastic waves in a fluid-saturated porous solid, I. Low-frequency range, J. Acoust. Soc. Am., 28, 168 – 178.
Biot, M. A. (1956b), Theory of propagation of elastic waves in a fluid-saturated porous solid, II. High-frequency range, J. Acoust. Soc. Am., 28, 179 – 191.
Central Geological Survey (1994), The final survey report of the Cho-Shui alluvial fan (in Chinese), Minist. of Econ. Affairs, Taipei, Taiwan. Central Geological Survey (1999), The hydrogeological survey report of the
Cho-Shui alluvial fan (in Chinese), Minist. of Econ. Affairs, Taipei, Taiwan.
Chia, Y. P., Y. S. Wang, H. P. Wu, C. J. Huang, C. W. Liu, M. L. Lin, and F. S. Jeng (2000), Changes of groundwater level in response to the Chi-Chi earthquake, in Proceedings of International Workshop on Annual Commemoration of Chi-Chi Earthquake, vol. I, Science Aspect, edited by C. H. Loh and W. I. Liao, pp. 317 – 328, Natl. Cent. for Res. on Earthquake Eng., Taipei, Taiwan.
Chia, Y. P., Y. S. Wang, J. J. Chiu, and C. W. Liu (2001), Changes of groundwater level due to the 1999 Chi-Chi earthquake in the Choshui River alluvial fan in Taiwan, Bull. Seismol. Soc. Am., 91, 1062 – 1068. Coleman, B. D., and W. Noll (1963), The thermodynamics of elastic
materials with heat conduction and viscosity, Arch. Rational Mech. Anal., 13, 167 – 178.
Department of Land Administration (2000), The records of the first and the second class GPS control stations, digital records, Minist. of the Interior, Taipei, Taiwan.
Drew, D. A., and S. L. Passman (1999), Theory of Multicomponent Fluids, 308 pp., Springer-Verlag, New York.
Eringen, A. C. (1980), Mechanics of Continua, 592 pp., R. E. Krieger, Huntington, N. Y.
Farrell, O. J., and B. Ross (1971), Solved Problems: Gamma and Beta Functions, Lendre Polynomials, Bessel Functions, 410 pp., Dover, Mineola, N. Y.
Flu¨gge, W. (1975), Viscoelasticity, 194 pp., Springer-Verlag, New York. Ge, S., and S. C. Stover (2000), Hydrodynamic response to strike- and
dip-slip faulting in half space, J. Geophys. Res., 105(B11), 25,513 – 25,524.
Grecksch, G., F. Roth, and H. J. Kumpel (1999), Cosesmic well level changes due to 1992 Roermond earthquake comparing to static deforma-tion of half space soludeforma-tions, Geophys. J. Int., 138, 470 – 478.
Green, A. E., and P. M. Naghdi (1965), A dynamical theorem of interacting continua, Int. J. Eng. Sci., 3, 231 – 241.
Hassanizadeh, M., and W. G. Gray (1979a), General conservation equations for multi-phase systems: 1. Averaging procedure, Adv. Water Resour., 2, 131 – 144.
Hassanizadeh, M., and W. G. Gray (1979b), General conservation equations for multi-phase systems: 2. Mass, momenta, energy, and entropy equa-tions, Adv. Water Resour., 2, 191 – 203.
Hassanizadeh, M., and W. G. Gray (1980), General conservation equations for multi-phase systems: 3. Constitutive theory for porous media flow, Adv. Water Resour., 3, 25 – 40.
Hsu, S. K. (1998), Plan for a groundwater monitoring network in Taiwan, Hydrogeol. Res., 6, 405 – 415.
Huang, B. S. (2000), Two dimensional reconstruction of the surface ground motion of an earthquake: The September 21, 1999, Chi-Chi, Taiwan earthquake, Geophys. Res. Lett., 27, 3025 – 3028.
Jeng, D. S., and T. L. Lee (2001), Dynamic response of porous seabed to ocean waves, Comput. Geotech., 28, 99 – 128.
Lamb, S. H. (1945), Hydrodynamics, 738 pp., Dover, Mineola, N. Y. Lee, M., T. K. Liu, K. F. Ma, and Y. M. Chang (2002), Coseismic
hydro-logical changes associated with dislocation of the September 21, 1999 Chichi earthquake, Taiwan, Geophys. Res. Lett., 29(17), 1824, doi:10.1029/2002GL015116.
Lin, Y. P., C. C. Lee, and Y. C. Tan (2000), Geostatistical approach for identification of transmissivity structure at Dulliu area in Taiwan, Environ. Geol., 40, 111 – 120.
Lin, Y. P., Y. C. Tan, and S. Rouhani (2001), Identifying spatial character-istics of transmissivity using simulated annealing and kriging methods, Environ. Geol., 41, 200 – 208.
Liu, S. Y., T.-C. J. Yeh, and R. Gardiner (2002), Effectiveness of hydraulic tomography: Numerical and sandbox experiments, Water Resour. Res., 38(4), 1034, doi:10.1029/2001WR000338.
Ma, K. F., C. T. Lee, and Y. B. Tsai (1999), The Chi-Chi, Taiwan earth-quake: Large surface displacement on an inland thrust fault, Eos Trans. AGU, 80, 605, 611.
Ma, K. F., J. Mori, S. J. Lee, and S. B. Yu (2001), Spatial and temporal distribution of slip for the 1999 Chi-Chi, Taiwan, earthquake, Bull. Seis-mol. Soc. Am., 91, 1069 – 1087.
Ohno, M., and H. Wakita (1997), A water well sensitive to seismic waves, Geophys. Res. Lett., 24, 691 – 694.
Prevost, J. H. (1982), Nonlinear transient phenomena in saturated porous media, Comput. Methods Appl. Mech. Eng., 20, 3 – 18.
Roeloffs, E. A. (1996), Poroelastic techniques in the study of earthquake-related hydrologic phenomena, Adv. Geophys., 37, 135 – 195.
Roeloffs, E. A. (1998), Persistent water level changes in a well near Park field, California, due to local and distant earthquakes, J. Geophys. Res., 103, 869 – 889.
Wang, C. C. (1970a), A new presentation theorem for isotropic functions, part I, Arch. Rational Mech. Anal., 36, 166 – 197.
Wang, C. C. (1970b), A new presentation theorem for isotropic functions, part II, Arch. Rational Mech. Anal., 36, 198 – 233.
Wang, C. Y., L. H. Cheng, C. V. Chin, and S. B. Yu (2001), Coseismic hydrologic response of an alluvial fan to the 1999 Chi-Chi earthquake, Taiwan, Geology, 29, 831 – 834.
Water Conservancy Agency (1999a), The annual report of the Groundwater Monitoring Network System (GMNS) records in Taiwan (in Chinese), Minist. of Econ. Affairs, Taipei, Taiwan.
Water Conservancy Agency (1999b), The annual report of the daily flow records in Taiwan (in Chinese), Minist. of Econ. Affairs, Taipei, Taiwan.
Water Conservancy Agency (2000), Analysis for the changes of surface levels and groundwater levels caused by the 921 Chi-Chi earthquake (in Chinese), Minist. of Econ. Affairs, Taipei, Taiwan.
Yeh, T.-C. J., and S. Y. Liu (2000), Hydraulic tomography: Development of a new aquifer test method, Water Resour. Res., 36, 2095 – 2105. Yu, S. B., et al. (2001), Preseismic deformation and coseismic displacement
associated with the Chi-Chi, Taiwan, earthquake, Bull. Seismol. Soc. Am., 91, 995 – 1012.
C.-H. Chen, Department of Civil Engineering, ChungKuo Institute of Technology, Taipei 106, Taiwan.
Y.-B. Lin, Council of Agriculture, Executive Yuan, Taipei 100, Taiwan. C.-W. Liu and Y.-C. Tan, Department of Bioenvironmental Systems Engineering, National Taiwan University, Taipei 10617, Taiwan. (yctan@ ccms.ntu.edu.tw)
T.-C. J. Yeh, Department of Hydrology and Water Resources, University of Arizona, Tucson, AZ 85721, USA.
