Visualization of aquifer response to periodic forcing

Advances in Water Resources 28 (2005) 819–834 www.elsevier.com/locate/advwatres

Visualization of aquifer response to periodic forcing Anthony J. Smith

a,*

, Lloyd R. Townley b, Michael G. Trefry

a

a b

CSIRO Land and Water, Private Bag No. 5, Wembley WA 6913, Australia Townley & Associates Pty Ltd., P.O. Box 3150, Nedlands WA 6909, Australia

Received 28 April 2004; received in revised form 19 January 2005; accepted 3 February 2005 Available online 9 April 2005

Abstract Many studies of periodic forcing and response in aquifers have focused on describing the induced fluctuations in hydraulic head, without much consideration of the time-dependent flows. Visualization techniques presented in this paper can be applied to obtain a more physically intuitive impression of groundwater motion in aquifers that undergo periodic fluctuations of hydraulic head. The concepts of velocity ellipse and displacement ellipse are introduced as methods for visualizing oscillatory flows associated with individual forcing modes. Cyclical trajectories illustrate the potential complexity of flow paths that can arise due to superposition of modal responses. The full periodic motion that results due to superposition of the mean flow and modal flows is visualized using streaklines. An animated time series of streaklines provides an intuitive impression of the flow and affords insight into apparent dispersion phenomenon that can arise due to periodic fluctuations in both the strength and direction of groundwater flow. Electronic animations are available from the authors.
2005 Elsevier Ltd. All rights reserved. Keywords: Visualization; Periodic; Cyclic; Aquifer; Groundwater; Streakline; Pathline

1. Introduction Periodic fluctuations of hydraulic head are common in many groundwater systems because periodic processes fundamentally drive the hydrological cycle. Daily and seasonal variations in the earthÕs radiation budget are reflected by periodic patterns of rainfall, evapotranspiration and groundwater recharge, which themselves drive periodic responses in aquifers [18,19,21,7,33]. Hydrologists have devoted considerable effort to studying the influence of surface water fluctuations on groundwater levels in contiguous aquifers, including effects from ocean tides and fluctuation of river stage

*

Corresponding author. Fax: +61 8 93336211. E-mail address: [email protected] (A.J. Smith).

0309-1708/$ - see front matter
2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2005.02.001

[11,14,26,6,23,33,35,34,2]. Quasi-cyclic variations in barometric pressure and fluctuation of external loads on aquifers due to ocean and earth tides also cause groundwater levels to oscillate [5,8,15,16]. Other fluctuations in groundwater response are induced by human activities such as diurnal and seasonal pumping from aquifers to meet irrigation and public water supply demands. Aquifer responses to periodic forcing are typically described in terms of the fluctuations that are induced in hydraulic head, without much consideration of the time-dependent flows. For example, in analysing twodimensional periodic flow fields, separate contour plots of the steady-state heads and the amplitudes and phases of head are commonly presented, e.g. [7,33,30]. Whereas this approach is relatively straightforward, it has limited utility because the time-dependent behavior of hydraulic head is not obvious from visual inspection of the plots. In addition, such plots do not provide a useful

820

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

impression of the flow. When more than one frequency is present, the interpretation of hydraulic head and flow becomes more difficult because two additional plots per frequency are required—one for amplitude and one for phase. These techniques are impractical for visualizing periodic flow in three-dimensional domains. Applying DarcyÕs Law to hydraulic heads that are periodic in time yields velocities that are also periodic. Visualizing such velocities by using separate plots of the steady-state component and of the amplitudes and phases of the periodic component can achieve a limited impression of the flow but the result is generally unsatisfying because it gives little insight about advective flowpaths and mixing. An alternate approach is to compute the distribution of total head at selected times throughout the period of fluctuation and then display these results as a time series of instantaneous equipotentials, e.g. [37]. By overlaying the flow vectors onto the equipotentials, a better impression of the groundwater motion can be achieved, but once again the result is inadequate for investigating advection and mixing. The purpose of this paper is to introduce simple visualization techniques that provide physically intuitive impressions of advective groundwater motion in periodic flow fields. Developed originally to visualize periodic fluctuations in the groundwater capture and release zones of lakes [31], these methods are generalized here to include three-dimensional flow and multiplefrequency forcing. The paper begins by considering the advective motions that arise due to spatially-uniform periodic flow and then extends this analysis to consider non-uniform periodic flow using a hypothetical groundwater flow problem. Animations for all examples presented in the paper, including supporting electronic materials, are available from the authors upon request.

odically in time. At each point in the flow field the velocity vector is described by n velocity components that each consists of a steady-state term and a cyclic term with m modes. This leads to n · (m + 1) terms to describe the total velocity at each point in the flow field. In this paper, attention is restricted to three-dimensional and two-dimensional flow fields. While recognising that spatially-uniform groundwater flow fields are physically impossible, they are considered, here, for theoretical interest and to introduce the visualization approaches. For three-dimensional flow m X uðtÞ ¼ us þ uc ðtÞ ¼ us þ uja cosðxj t  ujh Þ j¼1

vðtÞ ¼ vs þ vc ðtÞ ¼ vs þ

m X

vja cosðxj t  vjh Þ

ð1Þ

j¼1

wðtÞ ¼ ws þ wc ðtÞ ¼ ws þ

m X

wja cosðxj t  wjh Þ

j¼1

where u, v, and w are scalar Cartesian components of velocity [LT1] in x, y and z directions, respectively. Other variables are time t [T], angular frequency x = 2pf [T1], frequency f = 1/P [T1] and period P [T]. Subscript s denotes a steady-state term, c denotes a cyclic term, a denotes amplitude, h denotes phase, and superscript j denotes mode. For groundwater flow fields, u, v and w represent mean porewater velocities. Scalar components of displacement in the x, y and z directions are found by integrating Eq. (1) with respect to time over the arbitrary interval t ! t + dt m X uja X ðt; dtÞ ¼ us dt þ sinðxj dt þ xj t  ujh Þ xj j¼1 uja sinðxj t  ujh Þ xj m X vja sinðxj dt þ xj t  vjh Þ Y ðt; dtÞ ¼ vs dt þ j x j¼1 

2. Terminology

ð2Þ

The following terminology is adhered to throughout the text. A frequency component or mode is an individual sinusoid that is uniquely denned by its amplitude, period and phase. The term cyclic refers to purely oscillatory, time-dependent behaviour that can be expressed as a linear combination of modes. Steady-state describes the long-term average, or mean, value about which cyclical fluctuations occur. Together, the steady-state and cyclic terms add to give a periodic quantity. Other terms and symbols are introduced as they appear in the text.

These equations describe the total displacement of a fluid ÔparticleÕ, as measured from its starting position at time t to its position at time t + dt.

3. Advection in spatially-uniform periodic flow

4. Modal velocity ellipse

Consider an n-dimensional velocity field in which the independent velocity vector components fluctuate peri-

Each modal velocity vector has three scalar components, corresponding to the right-hand terms in

vj  aj sinðxj t  vjh Þ x m X wja sinðxj dt þ xj t  wjh Þ Zðt; dtÞ ¼ ws dt þ j x j¼1 

wja sinðxj t  wjh Þ xj

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

(a)

tion for the modal velocity ellipse can be written [40]

(b)

u - solid v - short dash w - long dash

Au2 þ Buv þ Cv2 ¼ 1

t

w

(0,0)

v u

(c)

Y X

Fig. 1. (a) Velocity components, (b) velocity ellipse and (c) displacement ellipses in three-dimensional spatially-uniform flow; ua = 1, va = 1, wa = 0.5, uh = 0, vh = p/2 and wh = p.

Eq. (1). Consider the single mode case m = 1. Then, neglecting the steady component for the present, we have uj ðtÞ ¼ uja cosðxj t  ujh Þ vj ðtÞ ¼ vja cosðxj t  vjh Þ

ð4Þ

with coefficients 1 A¼ j 2 2 ðua Þ sin ðsj Þ cosðsj Þ B¼ j j 2 ua va sin ðsj Þ 1 C¼ j 2 2 ðva Þ sin ðsj Þ

ð5Þ

where sj ¼ ujh  vjh is the phase difference between the component velocities. Eq. (4) describes an ellipse with principal axes rotated counter clockwise through the angle a, where a is given by   1 B 1 a ¼ tan ð6Þ 2 CA

Z

wj ðtÞ ¼ wja cosðxj t 

821

ð3Þ

wjh Þ

A parametric plot of these components (Fig. 1b) provides a simple graphical impression of the advective motion due to the modal response. The resultant ellipses are a type of Lissajous Figure or Bowditch Curve, which are well known in electronics [4] and the study of waves [25]. If more than one mode is present, i.e., m > 1, then each mode will be represented by a different velocity ellipse with its own size and orientation. A similar graphical representation of horizontal tidal currents is used in the study of ocean dynamics [12]. Whilst this approach has not previously been used by groundwater hydrologists, Meyboom [22] described a ‘‘gradient loop’’, which was obtained by plotting the horizontal groundwater gradient against the vertical gradient measured at the same location. The resultant plot was roughly elliptical and reflected distinct periodic forcing of groundwater flow due to seasonal water level fluctuations in an adjacent pond. Kacimov et al. [20] produced pathlines for particles advected in temporally-periodic velocity fields with single-mode, cyclostationary forcing. They demonstrated elliptical particle trajectories in response to sinusoidal forcing in homogeneous media, and the distortion of these pathlines in non-homogeneous media. Restricting attention to spatially-uniform, twodimensional flow in the x–y plane, the Cartesian equa-

Slopes and lengths of the major and minor axes can be found from the relations qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ðC  AÞ  ðA  CÞ þ B2 smaj ¼ qBffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð7Þ 2 ðC  AÞ þ ðA  CÞ þ B2 smin ¼ B sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ s2maj lmaj ¼ 2 A þ smaj ðB þ Csmaj Þ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð8Þ 1 þ s2min lmin ¼ 2 A þ smin ðB þ Csmin Þ Thus, modal velocity ellipses can be constructed either parametrically using Eq. (3) or geometrically by scaling and rotating a unit circle using the scaling factors lmaj and lmin and rotation angle a. Special attention is required in evaluating Eq. (5) through Eq. (8) due to discontinuities at sj = np, where n = 0, 1, 2, 3, . . ., 1. For the even-valued cases n = 0, 2, 4, 6, the length of the minor axis is zero and the ellipse reduces to a straight line, with motion occurring backward and forward along the direction of the major axis. In these cases, the slope and length of the major axis are found from sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 1 1 1 1 lim smaj ¼ ua va 2  2 þ þ 4 2 2 ; 4 sj !np=2 ua v a ua v a ua v a n ¼ 0; 4; 8; . . . lim smaj ¼ ua va

sj !np=2

n ¼ 2; 6; 10; . . . qffiffiffiffiffiffiffiffiffiffiffiffiffiffi lmaj ¼ 2 u2a þ v2a ;

1 1  2þ 2 ua v a

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi! 1 1 1 þ  ; u4a v4a u2a v2a

n ¼ 0; 2; 4; . . .

ð9Þ

822

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

For the remaining odd-valued cases, n = 1, 3, 5, . . ., the coefficients in Eq. (4) reduce to A ¼ 1=u2a , B = 0 and C ¼ 1=v2a , and Eq. (4) simplifies to the standard equation for a non-rotated ellipse with semimajor axis, ua, and semiminor axis, va u2 v 2 þ ¼1 u2a v2a

ð10Þ

The lengths of the major and minor axes are then lmaj = 2ua and lmin = 2va.

in the region of a particles motion have similar size and orientation. Irrespective of whether this condition is met, plotting displacement ellipses provides a useful technique for visualizing the spatial pattern of modal response. In spatially-uniform, two-dimensional flow in the x–y plane, the Cartesian equation for the modal displacement ellipse is Dx2 þ Exy þ Fy 2 ¼ 1

ð12Þ

with coefficients 5. Modal displacement ellipse

D¼

The modal components of displacement corresponding to the modal velocities are given by the right-hand terms in Eq. (2)

E¼ F ¼

uja uj sinðxj dt þ xj t  ujh Þ  aj sinðxj t  ujh Þ j x x j j v v j Y j ðtÞ ¼ aj sinðxj dt þ xj t  vh Þ  aj sinðxj t  vjh Þ x x j w wj Z j ðtÞ ¼ aj sinðxj dt þ xj t  wjh Þ  aj sinðxj t  wjh Þ x x X j ðtÞ ¼

ð11Þ Because there are now two parameters, t and dt, a parametric plot yields a family of displacement ellipses (Fig. 1c). These have the same shape and orientation as the modal velocity ellipse but are transformed in dimension by the scale factor 1/x [T]. Each displacement ellipse corresponds to the trajectory of an advected particle released at a different starting time t. All trajectories intersect at the point (X, Y, Z) = (0, 0, 0) because displacement for each mode is always equal to zero at travel times equal to integer multiples of the modal period, irrespective of the starting time. Displacement ellipses and modal trajectories are only coincident in spatially-uniform velocity fields, which are physically implausible in groundwater systems. As an advected particle moves from one point to another in a non-uniform velocity field, the modal velocity at its new location will be different to the modal velocity at its previous location. The true modal trajectory depends on the spatial distribution of velocity throughout the region that the particle travels, implying that particles have non-zero net displacement each period of travel—accumulating as particle drift. Shum [29] examined this phenomenon in porewater trajectories below a rippled water-sediment interface subject to cyclic wave action. Kacimov et al. [20] illustrated similar fluid advec-tion patterns in porous media with ‘‘parquettype’’ conductivity distributions and subject to cyclical forcing. Despite this limitation in non-uniform flow, a displacement ellipse provides a reasonable approximation of the true modal trajectory if all velocity ellipses

x2 ðuja Þ2 sin2 ðsj Þ x2 cosðsj Þ uja vja sin2 ðsj Þ x2

ð13Þ

ðvja Þ2 sin2 ðsj Þ

Note that D, E and F are simply the coefficients A, B and C from Eq. (5) with the amplitude terms scaled by 1/x. The angle of rotation and slopes and lengths of the principle axes can be calculated from Eqs. (6)–(10) by replacing the coefficients A, B and C by D, E and F, respectively.

6. Cyclical trajectory A fluid particleÕs cyclic trajectory is simply the superposition of its modal trajectories m X uja X c ðt; dtÞ ¼ sinðxj dt þ xj t  ujh Þ j x j¼1 uja sinðxj t  ujh Þ xj m X vja sinðxj dt þ xj t  vjh Þ Y c ðt; dtÞ ¼ j x j¼1 

ð14Þ

vj  aj sinðxj t  vjh Þ x m X wja sinðxj dt þ xj t  wjh Þ Z c ðt; dtÞ ¼ j x j¼1 

wja sinðxj t  wjh Þ xj

Parametric plots of the cyclic terms are potentially more complex than the component modal trajectories. Cyclic trajectories repeat the same pattern over either the longest modal period or possibly a longer period if a beat frequency arises due to superposition of modes with similar frequencies. This interference effect is akin to the audible ÔbeatingÕ that is heard when two sound waves of similar frequency are played simultaneously. In oceans, for example, tides with similar magnitude O1

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

(period = 1.076 days) and K1 (period = 0.997) constituents beat with an approximate 14-day period. Here, we refer to the period of repeat as the cyclic period and denoted it Pc, where xc = 2p/Pc. Fig. 2 depicts an example cyclic trajectory in spatially-uniform threedimensional flow with bi-modal forcing. The velocity scalar components are each composed of two modes with periods that differ by a factor 2 (e.g. semi-diurnal tide). The cyclic trajectory exhibits a double-loop pattern that repeats itself over the longer of the two modal periods, which is also the cyclic period in this example. If the period of one mode were three times the other, then the motion would be characterized by a triple-loop pattern, and so on. Potentially more complicated trajectories arise when the modal periods give rise to a beat frequency, as illustrated in Fig. 3. In this example, the modal periods differ by a factor 1.2 and the cyclic trajectory traces a more convoluted pattern that repeats itself over a cyclic period that is six times the smaller modal period and five times the larger modal period. In such cases the cyclic period is equal to the beat period given by Pc = P1P2/jP1  P2j.

(a)

(b)

u - solid v - short dash w - long dash

Y

t

X

Z

Fig. 2. (a) Velocity components and (b) cyclic trajectory in threedimensional spatially-uniform flow with bi-modal forcing; u1a = 3, v1a = 1.2, w1a = 1.5, u2a = 0.5, v2a = 0.5, w2a = 1, u1h = p/2, v1h = p, w1h = p, u2h = p, v2h = p/2, w2h = p/2.

(a)

(b)

u - solid v - short dash w - long dash

Y

t

X

823

7. Pathlines The full time-dependent motion of an advected fluid particle in spatially-uniform periodic flow is described by Eq. (2). Fig. 4 depicts pathlines for two-dimensional flow with single-mode forcing. In each example, the velocity ellipse is drawn centered on the point (us, vs), which is at the tip of the steady-state velocity vector. Three characteristic pathline shapes can be identified; waves, cusps and loops. In Example (a), the plot origin (0, 0) is inside of the velocity ellipse and the pathline is looped because both of the velocity components reverse direction at the same time, causing the pathline to double back upon itself. Example (b) depicts a cusp-shaped pathline that occurs only if the ellipse passes exactly through the plot origin. Cusps on the pathline indicate times at which the velocity is instantaneously zero, which correspond to the particle being at instantaneous stagnation points in the velocity field. In Example (c), the origin lies outside of the velocity ellipse and the generated pathline is smooth and wavy. Thus, the geometric relationship between the velocity ellipse and plot origin provides qualitative information about the type of advective motion. A special situation exists if the steady velocity is either zero or very small. If zero, the center of the velocity ellipse is located at the plot origin and the pathline traces the same shape as the displacement ellipse, with zero net displacement after each period of travel. If the steady velocity is small compared to the amplitude of the modal velocity, as in Example (d), then the pathline is characterised by large loops and fluid particles have small net displacements after each period of travel. In three-dimensional flow with multiple modes, the forms of loops, cusps and waves along a pathline are potentially more complicated. Figs. 5 and 6 depict the pathlines corresponding to the examples from Figs. 2 and 3 with the addition of steady-state velocity components. In Fig. 5, the double-loop pattern depicted in Fig. 2 is stretched out along the average flow path and the pathline repeats the same pattern over the cyclic period. The pathline in Fig. 6 is more intricate and exhibits three-dimensional forms of loops, cusps and waves along a single trajectory. Fig. 7 illustrates the potential complexity of pathlines in three-dimensional flow with bi-modal forcing.

Z

8. Streakline animation Fig. 3. (a) Velocity components and (b) cyclic trajectory in threedimensional spatially-uniform flow with bi-modal forcing and lowfrequency beat; u1a = 3, v1a = 1.2, w1a = 1.5, u2a = 0.5, v2a = 0.5, w2a = 1, u1h = 0, v1h = p/2, w1h = p, u2h = p/2, v2h = p/2, w2h = p/2.

Less familiar than the concept of a pathline is a streakline [36,3]. A streakline joins the instantaneous positions of all fluid particles that have passed through a single fixed location at some previous time. Consider

824

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

Velocity Components

Velocity Ellipse

Pathline

Y

v

u - solid; v- short dash

(a)

X

t

u- solid;

Y

v

v- short dash

(b)

u

X

t

u

v

u - solid; v- short dash

Y

(c)

X

t

u

Y

v

u - solid; v- short dash

(d)

t

X

u

Fig. 4. Periodic pathlines in two-dimensional spatially-uniform flow; (a) looped pathline, (b) cusp-shaped pathline, (c) wavy pathline and (d) oscillatory pathline.

(a)

(b)

(b)

(a) u - solid v - short dash w - long dash

u - solid v - short dash w - long dash

Y

t

Y

X

Z

t

X

Z

Fig. 5. (a) Velocity components and (b) pathlines in three-dimensional spatially-uniform flow with bi-modal forcing; us = 0.35, vs = 0.35, ws = 0.35, u1a = 3, v1a = 1.2, w1a = 1.5, u2a = 0.5, v2a = 0.5, w2a = 1, u1h = p/2, v1h = p, w1h = p, u2h = p, v2h = p/2, w2h = p.

Fig. 6. (a) Velocity components and (b) pathlines in three-dimensional spatially-uniform flow with bi-modal forcing and low-frequency beat; us = 0.35, vs = 0.35, ws = 0.35, u1a = 3, v1a = 1.2, w1a = 1.5, u2a = 0.5, v2a = 0.5, w2a = 1, u1h = 0, v1h = p/2, w1h = p, u2h = p/2, v2h = p/2, w2h = p/2.

the purely advective transport of a conservative tracer that is emitted from a continuous point source. The plume of tracer is reduced to a line, or streak, because there is no dispersion. Whereas a pathline represents the time history of a single particle, a streakline is generated at a particular instant in time and represents the

instantaneous positions of many fluid particles. Thus, pathlines and streaklines differ in transient flow. One approach to computing a streakline is to track a sequence of particles released at different starting times from a fixed point within the flow field. The positions of all the particles at a particular time of interest are

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

(a)

825

(b)

u - solid v - short dash w - long dash

X

t

Z

Y Fig. 7. (a) Velocity components and (b) pathlines in three-dimensional spatially-uniform flow with bi-modal forcing; us = 0.2, vs = 0.1, ws = 0.1, u1a = 1, v1a = 1, w1a = 3, u2a = 3, v2a = 1, w2a = 0.5, u1h = p, v1h = p/2, w1h = p/2, u2h = p/2, v2h = p, w2h = p.

joined together to form the streakline. If viewed in smooth succession, a time series of such streaklines produces an animation of the time-dependent flow. An advantage in analysing periodic flow is that streaklines need only be generated over one cyclic period since these can be animated in a loop to produce a continuous impression of the fluid motion. In spatiallyuniform periodic flow, streaklines can be computed analytically from the displacement equations (Eq. (2)). For example, to generate a streakline at time ts, particles are tracked from starting times t = ts  dt for a series of different travel times, say 0 6 t 6 tt, where tt is the time required for particles to travel completely through the aquifer from their point of release. The resultant particle locations at time ts constitute the streakline. Computation of streaklines in non-uniform periodic flow is handled in a similar way; however, the computational effort is far greater since a particle tracking code is required to compute particle displacements. Each point on a streakline is found as the endpoint of an individual pathline—terminated at time ts—and each streakline is generated as a series in time. If, for example, 100 points are required to give adequate resolution of a streaklineÕs geometry and the streakline is computed at 12 times within the cyclic period, then 1200 pathlines must be computed. If the flow field is visualized using 10 streaklines, then the number of pathlines to calculate is 12 000. To compute streaklines with suitable resolution, special care is required near regions of strongly divergent flow and in regions of high flows. Stretching in these areas can lead to a low density of points along a streakline and poor resolution of the streaklineÕs geometry. Remedying this requires an adaptive particle-tracking algorithm that dynamically allocates the travel times. Special handling of streaklines that are comprised of a number of separate pieces is also required in some instances. Fragmented streaklines occurs when some of the released particles end up outside of the flow domain and cannot be tracked further.

9. Example groundwater flow problem Groundwater velocity components that arise under field conditions conform to the physical laws that govern fluid flow through porous media. The above examples have considered advective motion in spatially-uniform, periodic velocity fields that were conceived without reference to physical problems. To demonstrate the visualization approaches in context of a non-uniform, physically-based flow field, a hypothetical groundwater flow problem with bi-modal forcing is now considered. The problem is intentionally contrived to focus attention on the visualization of the results rather than their practical significance. 9.1. Description The example problem (Fig. 8) consists of a rectangular flow domain (W · L) with spatially-uniform prescribed head boundary conditions along two opposing sides and streamlines along the remaining two sides. Steady, confined groundwater flow is driven by a head difference dh across the fixed head boundaries and periodic forcing is imposed by periodic injection and pumping. An injection well with injection period Pi is positioned mid-way along the left-hand streamline boundary and a pumping well with period Pp is located mid-way along the opposite streamline boundary. The steady injection and pumping rates are denoted Qis and Qps [L3T1] and the injection and pumping amplitudes Qia and Qpa [L3T1], respectively. The confined groundwater flow equation for twodimensional, transient flow through fully-saturated, inhomogeneous and anisotropic porous media is S

oh o2 h o2 h ¼ Tx 2 þ Ty 2 þ q ot ox oy

ð15Þ

where h is hydraulic head [L], T is aquifer transmissivity [L2T1], S is the bulk fluid-medium storage coefficient [1] and q is a general source term that includes pumping and injection [L/T]. The flow equation (15) assumes

826

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

Fig. 8. Geometry, boundary conditions and finite element mesh for the example groundwater flow problem.

that aquifer consolidation–expansion cycles are nonhysteretic, and that water is instantaneously released from and added to storage in response to periodic variation in head. The problem boundary conditions are hðx; 0; tÞ ¼ H hðx; L; tÞ ¼ H  dh qn ð0; y; tÞ ¼ 0; 0 6 y < L=2; L=2 < y 6 L qn ðW ; y; tÞ ¼ 0; 0 6 y < L=2; L=2 < y 6 L Qð0; L=2; tÞ ¼ Qis þ Qia cosðxi tÞ; y ¼ L=2

ð16Þ

QðW ; L=2; tÞ ¼ Qps þ Qpa cosðxp tÞ; y ¼ L=2 where qn is the component of groundwater flow normal to the boundary. Since the governing flow equation and boundary conditions are all linear, the superposition principle allows that the steady and modal responses can be solved independently and then recombined to obtain the full time-dependent solution. The example problem is fully described by twelve independent variables (L, W, T, S, ne, dh, Qis, Qps, Qia, Qpa, Pi, Pp) that are expressed in terms of two unit dimensions—length and time. This requires 10 (12  2) non-dimensional groups to describe the system behaviour, choices for these are listed and described in Table 1 along with the values adopted for the present example.

9.2. Computational approach Unlike wave phenomena, solutions to groundwater flow equations oscillate only when the forcing stresses are oscillatory. Wave propagation in a porous medium is normally assumed negligible because groundwater velocities and hence inertial forces are generally very small [9]. An aquifer that responds periodically to a periodic forcing behaves like an over damped forced oscillator [25,15,28]. Two approaches are used to analyse linear systems that experience forced periodic behaviour [17]. The first technique involves imposing a periodic forcing condition and solving an initial value problem that asymptotically approaches a periodic behavior—that is, a stationary state of constant amplitude and phase. The solution obtained consists of a periodic part and a transient part that dissipates with time [25]. In the second approach, the forcing conditions and system response are assumed to be fully periodic, implying that the transient term has vanished. A boundary value problem is solved for each forcing mode to obtain the head response in terms of amplitude and phase, e.g. [33,35]. The full time-dependent solution is then assembled by superposition of the single-mode solutions and steady solution. Bakker [2] recently applied the stationary state

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

827

Table 1 Non-dimensional ratios Group

Physical description

L W dh L Qis TW Qps TW Qia Qis

Aspect ration of flow domain

2

Regional hydraulic gradient

0.001

Qpa Qps L2 ne TP i

Value used

1

Ratio of the amplitude of pumping to the steady-state pumping rate

1

Non-dimensional response time that indicates the ability of the aquifer to propagate periodic forcing of period Pi—equal to the aquifer response time divided by the period of forcing Non-dimensional response time that indicates the ability of the aquifer to propagate periodic forcing of period Pp—equal to the aquifer response time divided by the period of forcing Aquifer elastic storage coefficient Effective porosity

L2 ne TP p S ne

approach to develop analytic element formulations of periodic, Dupuit–Forchheimer flow, which included periodic wells, line-sinks, and circular area-sinks. 9.3. Numerical solution using AQUIFEM-P The example problem in this paper was solved using the finite element modeling package AQUIFEM-P [32]. This package was designed specifically for simulating linear, periodically forced groundwater flow in twodimensional plan and section using the stationary state approach. Boundary forcing conditions are specified in the frequency domain and can include Type One (Dirichlet), Type Two (Neumann) and Type Three (Cauchy) boundaries, as well as distributed sources and sinks, point sources and sinks, and leakage from adjacent aquifers. The resultant steady (real valued) problem and modal (complex valued) problems are solved independently using the finite element method. Output from AQUIFEM-P includes (2 · m) + 1 nodal arrays consisting of the steady-state heads and the head amplitudes and head phases for each mode. Also outputs are (4 · m) + 1 element arrays consisting of the steady-state element velocities in each coordinate direction, and the velocity amplitudes and velocity phases for each mode in each coordinate direction. The full timedependent solution is obtained by superposition of the steady and modal responses; where, in general, a solution term fi at node or element i is evaluated as m X fi ðtÞ ¼ fsi þ faij cosðxj t  fhj Þ ð17Þ j¼1

0.0002

Non-dimensional transmissivity that indicates the ability of the aquifer to transport water at the steady-state injection rate Non-dimensional transmissivity that indicates the ability of the aquifer to transport water at the steady-state pumping rate Ratio of the amplitude of injection to the steady-state injection rate

0.0002

1.37

2.74

0.0001 0.25

Nodal mass balance terms are also output in the frequency domain, including storage fluxes, internal fluxes, boundary fluxes, source and sink fluxes, and leaky fluxes. A more detailed description of AQUIFEM-P and its use can be found in the ‘‘UserÕs Manual and Description’’ [32]. Smith [31] provides additional detail concerning the application of AQUIFEM-P to flow in vertical section. Validation of the code has included benchmarking against algebraic solutions for both one-dimensional and two-dimensional problems. Appendix A presents an example benchmark problem for the case of two-dimensional horizontal flow in a rectangular domain with boundaries constrained by single mode Dirichlet conditions. It is shown that AQUIFEM-P is quite accurate even for modest spatial discretizations. 9.4. Method for constructing ellipses To display the modal trajectories for the example problem, displacement ellipses for elements were constructed geometrically using Eqs. (6)–(10). Thus, the angle of rotation and lengths of the major and minor axes were computed from the velocity amplitudes and phases for each element. Ellipses were then drawn using the programming language PostScript [1] to scale and rotate a unit circle. 9.5. Method for particle tracking Although the steady and modal velocity fields computed by AQUIFEM-P are spatially discrete across the flow domain, the velocities within elements are

828

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

steady heads, (b and c) contours of head amplitudes for pumping and injection and (d and e) contours of head lags for pumping and injection. Lags in this instance are defined as the phase difference between forcing and response, divided by the period of forcing. Thus, head lags due to injection have values between zero and one and represent the fraction of the injection period that responses lag forcing at each location in the domain. Similarly, head lags due to pumping also have values between zero and one that are defined relative to the pumping period. From visual inspection of (b) and (c), the head responses to pumping and injection are both relatively small, with head amplitudes less than 0.2 at the injection and pumping wells. This implies that either the pumping and injection volumes are relatively minor components of the aquifer water balance or the aquifer response is manifested mostly as flows rather than water level fluctuations. In this example, it is known that the character-

spatially-uniform and continuous functions of time. Pathlines within elements were therefore computed analytically from the displacement Eq. (2). By knowing both the entrance point of a particle on the boundary of an element and the starting time for particle tracking within the element, the unique exit location from the element was straightforward to compute. Traversing the velocity field and searching for element entrance and exit locations was a convenient way to construct pathlines. The issue of velocity change during time steps, which arises if velocity is a function of discrete time, was avoided. The fact that exact, instantaneous locations of particles could be computed at any arbitrary time provided a distinct advantage for generating accurate streaklines. 9.6. Visualization of results Depicted in Fig. 9 are five typical plots of the steadystate and modal head distributions; (a) contours of (a)

(b)

(c)

0 .2 .04 .06 .08

.4

.04 .06 .08

.6

.02

.02

.8

1

5 6E-0

5E

(e)

5 -0 8E -04 1E -04 1.2E

1E-05 2E-05 3E-05

4E -05 -05

(d)

2E-05 4E-05 6E-05

Fig. 9. Periodic head distribution for the example groundwater flow problem. (a) Steady heads, (b) head amplitudes for injection, (c) head amplitudes for pumping, (d) head lags for injection and (e) head lags for pumping.

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

(a)

(b)

(c)

829

(d)

Fig. 10. Displacement ellipses and pathlines for the example groundwater flow problem. (a) Steady pathlines, (b) dispacement ellipses for injection, (c) displacement ellipses for pumping and (d) periodic pathlines.

istic aquifer response time, L2ne/T, [33] is relatively small compared to both the injection and pumping periods— that is, L2ne/TPi = 1.37 and L2ne/TPp = 2.74—which suggests that flow perturbations due to injection and pumping should propagate efficiently throughout the aquifer. In other words, the aquifer can respond almost in time with the forcing. From Fig. 9, it is clear that plots of head response do not provide an impression of these flows. Displacement ellipses and pathlines for the example problem are depicted in Fig. 10. Steady pathlines in (a) were generated by particle tracking forward in time from eight locations just inside of the up-gradient, prescribed head boundary. The same locations were used to generate periodic pathlines in (d) by releasing six particles from each location at equal time intervals within the cyclic period, which is Pi in this example. In (b) and (c), displacement ellipses for injection and pumping were drawn at gridded locations using the point-values of velocities. The ellipses are not true modal trajectories—because the flow is non-uniform—however they provide a useful impression of the modal response. Correct modal trajectories depend on the spatial distribution of velocity and have non-zero net displacement each period of travel. For example, Shum [29] illustrated this effect in porewater below a rippled water-sediment interface subject to periodic wave action. Displacement ellipses in Fig. 10 reveal a larger modal response to injection than pumping. This occurs because the forcing period for injection is twice the pumping period, allowing longer response to injection than pumping. In comparison, the plots of head amplitudes in Fig. 9 give the misleading impression that the aquifer responses to injection and pumping are almost identical. All displacement ellipses for both modes are closed, i.e., lmin = 0, indicating that the modal trajectories are mostly back and forth along straight lines; the conditions leading to this type of motion were discussed earlier in defining the velocity ellipse. Phase differences between the component velocities for each mode are not expected in this particular example because each

forcing mode is applied at only one location. Multiple forcing stresses of the same mode—for example, seasonal recharge and seasonal pumping—leads to open displacement ellipses such as those depicted in the earlier examples (e.g., Fig. 4). This effect is illustrated for the current example in Fig. 11, in which pumping and injection now have the same frequency but are out of phase by one quarter of a period. Open displacement ellipses with different shapes and orientations are evident due to the phase difference between forcing stresses and the induced phase lags within the aquifer. Note also that the responses to injection and pumping are now the same because Pi = Pp. Particle tracks in Figs. 10 and 11 illustrate that the steady pathlines are most perturbed by the modal responses where the direction of steady flow is orthogonal to the modal trajectories. Between the injection and pumping wells in Fig. 10, the steady and fluctuating flow directions are almost perpendicular, which leads to significant lateral perturbation of the ambient flow. Closer to the up-gradient and down-gradient prescribed head boundaries, the steady and fluctuating motions are almost aligned. The pathlines in these areas are not (a)

(b)

(c)

Fig. 11. Displacement ellipses and pathlines for a modified version of the example groundwater flow problem; pumping and injection have the same frequency but are one quarter of a period out of phase. (a) Steady pathlines, (b) displacement ellipses for injection and pumping and (c) periodic pathlines.

830

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

Fig. 12. Time sequence of streakline for the example groundwater flow problem.

wavy because fluctuations are back and forth along the mean flow direction. Whilst this type of pathline analysis is useful, it does not provide a complete impression of the timedependent flow. For example, it would be hard to ascertain the effects of pumping and injection on the shape and motion of a groundwater plume as it moved from up-gradient to down-gradient in the aquifer passed the injection and pumping wells. This difficulty arises because individual pathlines can cross due to different particles occupying the same location at different times. Pathlines fundamentally convey spatial information and therefore a different approach is required to visualize the flow dynamics. Fig. 12 presents a time sequence of streaklines that were generated at 12 times within the cyclic period. Each frame represents a snap shot of the flow taken at intervals 1/12th of the cyclic period. The streaklines can be likened to streaks of dye injected into the aquifer at

fixed-point locations and advected conservatively with the groundwater flow. Points of injection in Fig. 12 are indicated by dots. This static presentation of the streaklines can be enhanced by sequencing the frames to create an animation of the flow. The overall effect is to provide a dynamic impression of the release zone for injection and capture zone for pumping. For example, the fact that the capture zone of the pumping well undergoes two cycles of perturbation for each cycle of injection was not obvious from the pathline analysis. Furthermore, the two perturbation cycles in the capture zone of the pumping-well are slightly different from each other due to the effect of injection.

10. Summary and conclusions Periodic responses in aquifers are manifest both as fluctuations of hydraulic head and groundwater flow;

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

however, methods for visualizing periodic flows are not described in the groundwater literature. Because periodic flow responses are difficult to visualize based on conventional plots of the head response, simple methods for visualizing these flows are required. An improved understanding of aquifer response is achieved by exploring flow in addition to water levels. Plotting modal displacement ellipses is an effective approach for visualizing the modal responses. Conventional pathline analyses give a partial impression of the time-dependent flows but a much better impression is obtained from streakline plots that can be animated to provide a dynamic impression of the periodic groundwater motion. There is potential for applying the visualization methods described in this paper to explore apparent dispersion phenomenon in aquifers due to periodic fluctuations in the direction and strength of groundwater flow, e.g. [13,27,38,10,39]. Though not a theme of the paper, the examples presented here provide insight about the possible complexity of flow paths in periodic flow and the potential for apparent dispersion phenomenon in groundwater systems with distinct periodic forcing and response.

Acknowledgements This research was jointly supported by Murdoch University, Western Australia, and CSIRO Land and Water, Australia. The authors are grateful to Dr John H. Knight from CSIRO Land and Water for reviewing the draft manuscript. We also thank the three reviewers for AWR who provided valuable and constructive comments.

Appendix A. Validation of AQUIFEM-P against 2D analytical result In order to demonstrate the suitability of AQUIFEM-P in simulating stationary periodic flows, we first derive an analytical solution for a test problem, and then perform benchmarking runs to gauge accuracy of the numerical method. By performing a Fourier transformation of Eq. (15) a differential equation for the transformed head, U½hðx; y; tÞ   hðx; y; xÞ, is obtained: Dr2  hðx; y; xÞ ¼ ix hðx; y; xÞ

ðA:1Þ

where D = T/S is the aquifer diffusivity, T is the homogeneous aquifer transmissivity, S is the effective aquifer storativity, i is the imaginary unit, and U denotes the Fourier transformation operator. The groundwater head h(x, y, t) is recovered by inverse transformation of

831

hðx; y; tÞ ¼ U1 ½hðx; y; xÞ . Eq. (A.1) is a Helmholtz equation with a complex coefficient (ix). Consider the following problem. On the boundaries of the twodimensional domain [x, y] 2 [0, a] · [0, b], the solution hðx; y; xÞ is restricted by the Dirichlet boundary conditions hð0; y; xÞ ¼ h1 ðxÞ hða; y; xÞ ¼ h2 ðxÞ hðx; 0; xÞ ¼ h3 ðxÞ

ðA:2Þ

hðx; b; xÞ ¼ h4 ðxÞ h3 ðxÞ; The boundary condition functions h1 ðxÞ; h2 ðxÞ;  h4 ðxÞ are complex valued, with complex magnitudes corresponding to the amplitudes of the harmonic head conditions and complex arguments corresponding to the phases of the head conditions. The solution of this Helmholtz problem takes two parts. The first part is the formal solution of the sub-problem where three of the four boundary forcing functions in Eq. (A.2) are zero and the remaining boundary function is non-zero. This corresponds to a physical groundwater flow problem where three of the four fixed-head boundaries have head values that are constant in time, and the remaining boundary has a time-dependent fluctuation. The second part of the solution makes use of the linearity of the Helmholtz equation to construct by superposition a general solution to the case where all boundaries have nonzero forcing functions. Firstly, we take h1 ðxÞ ¼ h2 ðxÞ ¼ h3 ðxÞ ¼ 0 and h4 ðxÞ 6¼ 0. Writing hðx; y; xÞ ¼ X ðx; xÞY ðy; xÞ and substituting into Eq. (A.1) yields two separated ordinary differential equations: X 00 ðx; xÞ ¼ kX ðx; xÞ   ix Y 00 ðy; xÞ ¼    k Y ðy; xÞ D

ðA:3Þ

where the separation constant, k, is as yet unspecified. The fundamental solutions of the two equations are pffiffiffi pffiffiffi X ðx; xÞ ¼ A cosð kxÞ þ B sinð kxÞ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! ix ix Y ðy; xÞ ¼ E cos   ky þ F sin   ky D D ðA:4Þ The integration constants A, B, E, F and the separation constant k are to be fixed by application of the boundary conditions. By assumption, X(0; x) = X(a; x) = 0 which, pffiffiffi for non-trivial solutions, implies that A = 0 and k ¼ np=a for n = 1, 2, 3, . . . Similarly, Y(0; x) = 0 leads to E = 0. It is seen that infinitely many values of the k separation constant satisfy the boundary conditions. To

832

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

proceed we express the solution as a sum over these terms:  hðx; y; xÞ ¼

1 X n¼1

! rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

npx ix np2 y C n sin sin   a D a ðA:5Þ

where Cn is the (lumped) expansion coefficient of the nth term. The Cn coefficients may be fixed by applying the remaining boundary condition: hðx; b; xÞ ¼

1 X

C n sin

npx

n¼1

a

! rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ix np2 b sin   D a

¼ h4 ðxÞ

ðA:6Þ

Using the orthogonality of the sinðnpx=aÞ functions over x 2 [0, a], multiplying both sides of Eq. (A.6) by sinðnp=aÞ, integrating over the domain of x and rearranging finally yields 1  hðx; y; xÞ X ¼  h4 ðxÞ n¼1

n

2ð1  ð1Þ Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 b  np np sin  ix D a

 sin

npx a

! rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ix np2 y sin   D a

ðA:7Þ

This sum only contains odd terms since the coefficient numerator vanishes for even values of n. It is also clear that the solution scales with the boundary condition value  h4 ðxÞ. We note that corresponding solutions for cases with the only non-zero boundary condition on other boundaries can be gained either by reflecting the y coordinate, or swapping x with y and a with b, or both. The linearity of Eq. (A.1) ensures that the general solution to the multiple Dirichlet boundary problem is a superposition of these transformed solutions. The final result is then

hðx; y; x; a; b; D; h1 ; h2 ; h3 ; h4 Þ n 1 X 2ð1  ð1Þ Þ sin npy b   ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi np 2 ix n¼1 np sin D b a " ! rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

np2 ix  h1 sin ða  xÞ   D b !# rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ix np2  þ h2 sin x   D b n 1 X 2ð1  ð1Þ Þ sin npx a  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ

2 n¼1 np sin b  ix  np D b " ! rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

np2 ix  h3 sin ðb  yÞ   D a !# rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

np2 ix y þ h4 sin   D a

ðA:8Þ

Returning to the groundwater head, if the frequency spectrum of the head boundary conditions contains a set of discrete modes xm, then the solution for the head fluctuation in the domain is simply hðx; y; tÞ ¼ R

X

hðx; y; xm ; a; b; D; h1 ðxm Þ; h2 ðxm Þ;

m

h3 ðxm Þ; h4 ðxm ÞÞe

! ixm t

ðA:9Þ

A.1. Comparison with AQUIFEM-P The output of AQUIFEM-P was compared with the analytical result of Eq. (A.9) for the test case a = 1, b = 1.62, D = 2.62, x = 2p, h1 ¼ 0:025e2pi=5 ,  h2 ¼ 0:08epi=2 , h3 ¼ 0:1, h4 ¼ 0:05epi . Here the four boundaries are forced at the same frequency but with different amplitudes and phases. A triangular finite element grid

Fig. A.1. AQUIFEM-P grid and solution for the benchmark problem.

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834

Fig. A.2. Convergence of the percentage residual measures with respect to analytical series truncation limit.

with 3871 nodes was used to represent the problem domain. Comparison was made for a randomly chosen subset of kmax = 362 nodes within the domain. Boundary nodes were excluded since the analytical result uses a Fourier approach which is slowly convergent at the boundaries. Fig. A.1 shows the finite element grid used in the AQUIFEM-P simulation and the amplitude and phase maps of the analytical solution (comparison nodes highlighted). Agreement was measured by calculating the percentage normalized root-mean-square residuals for the amplitudes and means, i.e. " # 2 1=2 1 X ðjhnumerical j  jhanalytical jÞ Ramplitude ¼ 100 k max k jhnumerical j " #1=2 1 X ðargðhnumerical Þ  argðhanalytical ÞÞ2 Rphase ¼ 100 k max k argðhnumerical Þ ðA:10Þ where k indexes the comparison nodes. Fig. A.2 shows how the percentage residuals decline as the series truncation limits in Eq. (A.8) increase. The overall agreement between the AQUIFEM-P and analytical results is encouraging, especially for the amplitudes.

References [1] Adobe systems incorporated. PostScript language: tutorial and cookbook. Addison-Wesley: 1991. [2] Bakker M. Transient analytic elements for periodic Dupuit– Forchheimer flow. Adv Water Resour 2004;27:3–12. [3] Bear J. Dynamics of fluids in porous media. New York: Dover Publications; 1972. [4] Borowski EJ, Borwein JM. Collins reference dictionary of mathematics. London: Collins; 1989. [5] Bredehoeft JD. Response of well-aquifer systems to earth tides. J Geophys Res 1967;72(12):3075–87. [6] Carr PA, Van Der Kamp GS. Determining aquifer characteristics by the tidal method. Water Resour Res 1969;5(5):1023–31.

833

[7] Colville JS, Holmes JW. Water table fluctuations under forest and pasture in a Karstic region of Southern Australia. J Hydrol 1972;17:61–80. [8] Cooper HH, Bredehoeft JD, Papadopulos IS, Bennett RR. The response of well-aquifer systems to seismic waves. J Geophys Res 1965;70(16):3915–26. [9] Davis SN, De Wiest RJM. Hydrogeology. New York: John Wiley & Sons; 1966. [10] Farrell DA, Woodbury AD, Sudicky EA, Rivett MO. Stochastic and deterministic analysis of dispersion in unsteady flow at the Borden Tracer-Test site, Ontario, Canada. J Contamin Hydrol 1994;15:159–85. [11] Ferris JG. Cyclic fluctuations of water level as a basis for determining aquifer transmissibility. Internal Assoc Sci Hydrol 1952:148–55. General Assembly Brussels. [12] Foreman MGG, Henry RF. The harmonic analysis of tidal model time series. Adv Water Resour 1989;12:109–20. [13] Goode DJ, Konikow LF. Apparent dispersion in transient groundwater flow. Water Resour Res 1990;26(10):2339–51. [14] Gregg DO. An analysis of ground-water fluctuations caused by ocean tides in Glynn County. Georgia, Ground Water 1966;4(3): 24–32. [15] Hsieh PA, Bredehoeft JD, Farr JM. Determination of aquifer transmissivity from earth tide analysis. Water Resour Res 1987; 23(10):1824–32. [16] Hsieh PA, Bredehoeft JD, Rojstaczer A. Response of well aquifer systems to earth tides: problem revisited. Water Resour Res 1988; 24(3):468–72. [17] Hudson JA. The excitation and propagation of elastic waves. London: Cambridge University Press; 1980. [18] Jacob CE. Correlation of ground-water levels and precipitation on Long Island, New York: Part I—Theory. Trans, Am Geophys Union 1943:564–72. [19] Jacob CE. Correlation of ground-water levels and precipitation on Long Island, New York: Part II—Correlation of data. Trans Am Geophys Union 1944:928–39. [20] Kacimov AR, Obnosov YuV, Yakimov ND. Groundwater flow in a medium with a parquet-type conductivity distribution. J Hydrol 1999;226:242–9. [21] Maasland M. Water table fluctuations induced by intermittent recharge. J Geophys Res 1959;64(5):549–59. [22] Meyboom P. Unsteady groundwater flow near a willow ring in hummocky moraine. J Hydrol 1966;4:38–62. [23] Nevulis RH, Davis DR, Sorooshian S. Analysis of natural groundwater level variation for hydrogeologic conceptualisation, Hanford site, Washington. Water Resour Res 1989;25(7): 1519–29. [25] Pain HP. The physics of vibrations and waves. London: John Wiley & Sons; 1968. [26] Pinder GF, Bredehoeft JD, Cooper HH. Determination of aquifer diffusivity from aquifer response to fluctuations in river stage. Water Resour Res 1969;5(4):850–5. [27] Rehfeldt KR, Gelhar LW. Stochastic analysis of dispersion in unsteady flow in heterogeneous aquifers. Water Resour Res 1992; 28(8):2085–99. [28] Sears WS, Zemansky MW, Young HD. University physics. Massachusetts: Addison-Wesley Publishing Company; 1987. [29] Shum KT. Wave-induced advective transport below a rippled water-sediment interface. J Geophys Reas 1992;97(C1):789–808. [30] Sun H. A two-dimensional analytical solution of groundwater response to tidal loading in an estuary. Water Resour Res 1997;33(6):1429–35. [31] Smith AJ. Periodic forcing of surface water-groundwater interaction: modelling in vertical section. PhD, Murdoch University, Western Australia, 1999. [32] Townley LR. AQUIFEM-P: a periodic finite element aquifer flow model: userÕs manual and description (Version 1.0). CSIRO

834

[33] [34]

[35] [36]

A.J. Smith et al. / Advances in Water Resources 28 (2005) 819–834 Institute of Natural Resources and Environment, Division of Water Resources Technical Memorandum 93/13, 1993. Townley LR. The response of aquifers to periodic forcing. Adv Water Res 1995;18(3):125–46. Trefry MG, Bekele E. Structural characterization of an island aquifer via tidal methods, Water Resour Res 40(1):8 January 2004; W01505, 10.1029/2003WR002003. Trefry MG. Periodic forcing in composite aquifers. Adv Water Resour 1999;22(6):645–56. Vallentine HR. Applied hydrodynamics. 2nd ed. London: Butterworths; 1967.

[37] Vandenberg A. Regional groundwater motion in response to an oscillating water table. J Hydrol 1980;47:333–48. [38] Webster IT, Taylor JH. Rotational dispersion in porous media due to fluctuating flows. Water Resour Res 1992;28: 109–19. [39] Zhang D, Neuman P. Head and velocity covariances under quasisteady-state flow and their effects on advective transport. Water Resour Res 1996;32(1):77–83. [40] Zwikker C. The advanced geometry of plane curves and their applications, formerly titled: advance plane geometry. New York: Dover Publications; 1966.