Transient flow between aquifers and surface water : analytically derived field-scale hydraulic heads and fluxes

The increasing importance of catchment-scale and basin-scale models of the hydrological cycle makes it desirable to have a simple, yet physically realistic model for lateral subsurface water flow. As a first building block towards such a model, analytical solutions are presented for horizontal groundwater flow to surface waters held at prescribed water levels for aquifers with parallel and radial flow. The solutions are valid for a wide array of initial and boundary conditions and additions or withdrawals of water, and can handle discharge into as well as lateral infiltration from the surface water. Expressions for the average hydraulic head, the flux to or from the surface water, and the aquifer-scale hydraulic conductivity are developed to provide output at the scale of the modelled system rather than just point-scale values. The upscaled conductivity is time-variant. It does not depend on the magnitude of the flux but is determined by medium properties as well as the external forcings that drive the flow. For the systems studied, with lateral travel distances not exceeding 10 m, the circular aquifers respond very differently from the infinite-strip aquifers. The modelled fluxes are sensitive to the magnitude of the storage coefficient. For phreatic aquifers a value of 0.2 is argued to be representative, but considerable variations are likely. The effect of varying distributions over the day of recharge damps out rapidly; a soil water model that can provide accurate daily totals is preferable over a less accurate model hat correctly estimates the timing of recharge peaks.


Introduction
Many agriculturally productive regions in temperate climate zones are located in areas with little topography and shallow groundwater, such as delta areas.The precipitation surplus is often discharged via dense, partially man-made drainage systems (Lennartz et al., 2009).For individual fields, drainage theory based on analytical solutions for the predominantly horizontal flow in the phreatic aquifer has proven its value for several decades (Hooghoudt, 1940;Dumm, 1954;Kraijenhoff van de Leur, 1958; van Schilfgaarde, 1970).
Currently, the changing climate drives efforts to model the terrestrial hydrological cycle at the scale of entire catchments and basins.The horizontal saturated flow to and from the drainage network is an important segment of the hydrological cycle at these larger scales.Basin-scale and global models tend to focus on the vertical column covering the unsaturated zone and the atmosphere.In many models, lateral flows between columns are represented in a conceptual manner (see the overview by Nijssen et al., 2001;Samaniego et al., 2010, for a recent example).Often the focus is more on mountainous and hilly areas than on flatter terrain (see Gong et al., 2011).The model columns can have horizontal dimensions that render point-values of typical hydraulic parameters and variables such as hydraulic conductivity, hydraulic head, and flux density useless.Therefore, an approach is desirable that is based on a simplified flow description that better reflects the essential features of lateral subsurface flow than the conceptual approaches used so far and still expresses the results in terms of large-scale variables: the flux between the groundwater and the surface water, and an average measure of the phreatic level.
Published by Copernicus Publications on behalf of the European Geosciences Union.

G. H. de Rooij: Transient flow between aquifers and surface water
An alternative path is offered by highly advanced integrative modelling that couples processes at the soil surface, the soil, the groundwater, and the surface water (e.g.Hy-droGeoSphere) (Brunner and Simmons, 2011;Hydrogeosphere, 2012).Such models hold great potential for complicated local to continental studies (which will often involve solute transport) in order to improve management strategies or guide measures to protect groundwater and surface water (e.g.Li et al., 2008).Their sophistication makes them very data-intensive.The optimal use of these models requires a prolonged and dedicated effort to set-up the model, provide the input, and store and analyze its output.The coupling with atmospheric models is less advanced than for the large-scale models discussed above.In this manuscript we opt instead for a more explorative, less data-intensive, and computationally light approach.
While traditional drainage theory is of limited use for catchment-scale models, its analytical approach may support the development of a less conceptual and more physical representation of the fluxes between the groundwater and the surface water.De Rooij (2009) explored the upscaled equivalents of conventional Darcian flow descriptors by their energy-conserving volume averages.De Rooij (2011) recently showed that aquifer-scale steady-state horizontal saturated flows behave in a Darcian way in that the flux between the groundwater and the surface water is directly proportional to the difference between the energy-conserving averaged hydraulic heads of the two water bodies.Thus, in principle, the average groundwater level and the surface water level, together with an upscaled hydraulic conductivity would suffice to model groundwater-generated stream discharge.Such average groundwater levels and upscaled conductivities can readily be derived from the analytical solutions to the saturated flow problems treated in drainage theory.
Obviously, steady-state solutions will not be adequate for many practical problems.Based on de Rooij's (2011) proof of principle, this paper therefore explores linearized transient groundwater flows in order to examine parallel and radial flows toward or from surface waters.It does so through analytical solutions of the differential equations describing the flow.The advantages over numerical solutions are that the resulting expressions provide a more profound insight into the fundamental behaviour of the systems and that upscaled parameters and variables can be calculated exactly.This does not imply that future applications should necessarily be analytical also, but the insight gained from the analytically derived relationships can inform future implementations, be they analytical or numerical.
Generic solutions are developed that cover nine different scenarios that reflect combinations of different forcing mechanisms and changes in these forcings, caused, for instance, by the commencement and cessation of rainfall, or human manipulation of surface water levels.The term forcing in this paper refers to initial and boundary conditions, recharge rate, and head-dependent recharge.The solutions do not appear to have been published before.These solutions are used to analyse the behaviour of the aquifers under different conditions, and to compare the effects of parallel and radial geometry on the hydraulic head and the flow.Also, since precise rainfall predictions at the field scale are impossible, the effect of the temporal distribution on recharge (generated by infiltrating rainfall) is considered.
In view of potential applications in large-scale models that cannot accommodate local (point-scale) values of heads and fluxes, expressions are developed for the average hydraulic head and the flux at the groundwater-surface water interface.The relationship between the two is investigated in some detail.The solutions show that the linear relationship between average hydraulic head and steady-state discharge proved by de Rooij (2011) does not exist for transient flows.Instead, a more complicated, time-dependent, but still explicit relationship connects the two.This relationship allows the calculation of the hydraulic conductivity at the field scale (the scale of the system between the zero-flux boundary at the axis of symmetry and the surface water) expressed solely in terms of the initial and boundary conditions and the geohydrological properties of the subsurface.Furthermore, fluxes towards the surface water and average hydraulic heads can be calculated directly from the forcings and the geohydrological parameters.The theory developed here thus provides the building blocks for an approach that can connect predominantly horizontal, field-scale groundwater flows to the essentially vertical hydrology of soil-vegetation-atmosphere exchange processes.

Parallel flow: governing partial differential equation (PDE) and initial (IC) and boundary conditions (BC)
Invoking the Dupuit assumptions for groundwater flow eliminates vertical gradients in the hydraulic head, and only the horizontal coordinates remain.Phreatic aquifers may receive recharge from the unsaturated zone above that is independent of the local hydraulic head and may exchange water with a deeper aquifer if the separating aquitard is somewhat permeable.Such exchange fluxes are assumed here to be proportional to the local hydraulic head.For a uniform, isotropic, phreatic aquifer overlying a level aquitard, the governing PDE then becomes: where x 1 and x 2 [L] are the horizontal coordinates, t * is time [T], H is the hydraulic head [L], defined with respect to the top of the underlying aquitard, K [LT −1 ] is the hydraulic conductivity, R [LT −1 ] is the recharge or loss to evapotranspiration (R may be time dependent), µ is the storage coefficient (occasionally termed drainable porosity for a phreatic aquifer, e.g.van Schilfgaarde, 1974), and a [T −1 ] (≤0) and b [LT −1 ] are constants determining the exchange with the deeper aquifer (see also Appendix D for a list of symbols).If the deeper aquifer has a constant and uniform hydraulic head H 2 , −a −1 is the resistance of the aquitard, and b = −aH 2 .Equation ( 1) is the Boussinesq equation with additional production terms and does not have a general analytical solution.
To make the equation analytically tractable it needs to be linearized by assuming that the variation in H is small with respect to H , and that µ is a constant: where D [L] is the constant water level above the aquitard.
Figure 1 gives a definition sketch of the original and the linearized problem.The combination of the Dupuit assumptions and the linearization has a sound footing in classical drainage theory (e.g.Dumm, 1954;Kraijenhoff van de Leur, 1958; van Schilfgaarde, 1970;Wesseling, 1979).Van Schilfgaarde (1974) gives a thoughtful discussion of the assumptions underlying the above linearization.See Appendix A for a quantitative treatment of the storage coefficient.
To analyze flow towards parallel drains, ditches, or streams with spacing 2L [L], we drop the second horizontal coordinate since the flow lines are all perpendicular to it.We also make the independent variables dimensionless by the following transformations: to obtain: This equation needs to be solved for various cases, all in the domain 0 ≤ x < 1 and t > 0. For all cases, x = 0 lies at the midpoint between two surface water bodies or drains.It therefore constitutes an axis of symmetry where there is no flow: The BC at x = 1, the IC, and the values of a, b, and R, vary from case to case: Case 1.Initial hydrostatic equilibrium with a step change of H (1,0) at t = 0.This case reflects the sudden increase or decrease of the ditch water level, for instance to increase the groundwater level during dry periods.
IC, BC, and parameter values: Case 2. Initial hydrostatic equilibrium with step change of the water level in the soil at t = 0, reflecting a pulsed water input (e.g. by short, heavy rainfall).The surface water level remains constant.Mathematically, this problem is identical to Case 1, and the same solution applies.Case 3. Like Case 1, but with constant recharge or loss.This case reflects the sudden increase or decrease of the ditch water level, while a steady flux to or from the unsaturated zone to the phreatic aquifer is maintained.
IC, BC, and parameter values: Case 4. Like Case 2, but with constant recharge or loss.This reflects a pulsed water input followed by gentler recharge or loss.Mathematically, Case 4 is identical to Case 3. Case 5. Like Case 1, but with recharge or loss proportional (but not necessarily directly proportional) to the hydraulic head.This case reflects the sudden increase or decrease of the ditch water level, with a flux to or from the phreatic aquifer that consists of a constant and a H -dependent component.This flux can be composed of recharge from or flow to the unsaturated zone and to a deeper aquifer across an aquitard.
IC, BC, and parameter values: Case 6.Like Case 2, but with recharge or loss a linear function of H , like Case 5. Mathematically, this problem is identical to that of Case 5. Case 7. Initial hydrostatic equilibrium.Constant recharge (or loss) R 1 [LT −1 ] for 0 < t < t 1 , zero loss or recharge for t ≥ t 1 .This case can represent a single prolonged rainfall event.
IC, BC, and parameter values: Case 8. Like Case 7, but with constant recharge or loss R 2 [LT −1 ] for t ≥ t 1 .This can represent rainfall followed by constant (possibly potential) evapotranspiration.IC, BC, and parameter values: Case 9. Like Case 8, but with recharge/loss linearly varying with H , and with a sudden ditch water level change at t = 0.This case arises when the replenished phreatic aquifer loses water to the aquifer below, or when the delivery of water to the unsaturated zone is limited by the dropping hydraulic head in the phreatic aquifer.The BC at x = 1 adds additional flexibility, compared to Cases 7 and 8. IC, BC, and parameter values: Cases 1 through 6 are all covered by Eqs. ( 13)-(15) of Case 5, which is itself a special case of Case 9. Cases 7 and 8 are special cases of Case 9 as well, governed by Eqs. ( 24) through ( 27).Thus, all problems listed above are special cases of the solution to this most general case.The derivation of the solution for this case is presented in Appendix B. The resulting expression for the hydraulic head reads: where n is a counter, and u(t) is the Heaviside step function.

Radial flow: governing partial differential equation (PDE) and initial (IC) and boundary conditions (BC)
In axisymmetrical flows in a circular aquifer, the Dupuit assumptions are invoked again and both head-dependent andindependent exchanges of water with a deeper aquifer and/or the overlying unsaturated zone are permitted.For a uniform porous medium and IC and BCs that are independent of the location on the boundary, all flows will be radial, and the angular coordinate can be eliminated.The governing PDE then becomes: where r * [L] is the horizontal, radial coordinate.Linearizing as before by assuming the vertical extent of the saturated zone as well as µ constant gives: We seek solutions of this equation for finite radial domains (e.g., circular fields or reclaimed areas (polders) surrounded by a ditch): 0 ≤ r * < L. (The notation L is retained to facilitate comparison with the parallel flow solution.)Introducing dimensionless variables results in: (Non-consecutively numbered equations have been introduced before and are labelled by their original number).We can develop the same nine cases as for parallel flow, this time for the domain 0 ≤ r < 1 and t > 0. For all cases, r = 0 constitutes an axis of symmetry where there is no flow: The BC at r = 1, the IC, and the values of a, b, and R, vary from case to case.A very general problem analogous to Case 9 for parallel flow is defined by: Appendix C gives the derivation of the solution, which is: where J i (y) is the Bessel function of the first kind and order i, and α n are the roots of 2.3 Relationship between the flux to/from the surface water and the average hydraulic head

Parallel flow
In Eq. ( 28), only the cosine term depends on x.By expressing the gradient of the hydraulic head in dimensional form, it can be used to find the horizontal flux Q(x,t) [L 2 T −1 ] in the case of parallel flow: where x 2 [L] represents the horizontal coordinate running parallel to the aquifer boundary with the surface water.At the interface with the surface water (x = 1), this simplifies to The average hydraulic head H (t)in the aquifer is Note that for t = 0, the term in braces equals H 0 -H A and can be brought outside of the sum.Since we have for the series  (Euler, in Berggren et al., 2004, p. 116), the average head at t = 0 equals the initial head H 0 , which is correct.From Eqs. ( 38) and ( 39), the flux across a stretch x 2 [L] of the aquifer boundary is related to the average hydraulic head in the aquifer as: where For t = 0, all M n (0) are equal to H 0 -H A , and Q(1,0) is infinitely large when H 0 = H A and zero when H 0 = H A .
The non-linearity in the Q−H relationship of Eq. ( 41) arises from the term with the series.If only the first terms of both series are retained, the relationship becomes linear: The series does not converge fast for all t, however, so this simplification should be used with care.For t approaching infinity, the effect of the initial condition and of R 1 damps out, and only the BCs and R 2 affect the head.For a non-leaky aquifer without H -dependent evapotranspiration (a = b = 0), Eq. ( 41) for infinitely large time becomes where the values of the series were taken from Euler (in Berggren et al., 2004, p. 116).This result is consistent with de Rooij's (2011) steady-state analysis.Note that Eqs. ( 43) and ( 44) differ by 17 %.

Radial flow
In Eq. ( 35), only term with the Bessel functions describes the dependency of H on r.The horizontal flux Q(r,t) = −2πrKDL ∂H ∂r * [L 3 T −1 ] for radial flow follows from the gradient of the hydraulic head in dimensional form.To find this gradient, the following relationship is used (Abramowitz and Stegun, 1965, 9.1.30): From Eq. ( 35) then follows At the interface with the open water at r = 1 (r * = L) this gives If t = 0, Eq. ( 47) reduces to which equals infinity for et al., 1993, Eq. 2.7).The average hydraulic head is 16, 649-669, 2012 www.hydrol-earth-syst-sci.net/16/649/2012/ For t = 0, and a = b = R 1 = R 2 = 0, Eq. ( 50) reduces to With α 0 = 2.4048255577 (Abramowitz and Stegun, 1965, p. 409), the contribution of the first term is 0.1729 (69 % of the sum).With Eqs. ( 47) and ( 50), the Q−H relationship is where As in the case of parallel flow, the Q−H relationship can be linearized by retaining only the first term of both series in Eq. ( 52), but again, this approximation can be poor owing to slow convergence of the series.Note, incidentally, that the linearization of Eqs. ( 43) and ( 53) makes Q proportional to H − H A , which is the definition of a linear reservoir (e.g.Fenicia et al., 2006).The reservoir coefficients are fully defined in terms of porous media properties, aquifer dimensions, and aquifer geometry.

The upscaled hydraulic conductivity
The equation pairs ( 41)-( 42) and ( 52)-( 53) relate the flux across the groundwater-surface water interface to the difference between the average hydraulic heads on either side of this interface (in the surface water body, the average hydraulic head can be assumed identical to any local value).The proportionality constants K up [LT −1 ] and 2πLK ur [L 2 T −1 ] for parallel and radial flow are defined from Eqs. ( 41) and ( 52), respectively, as: The quantities K up and K ur represent field-scale equivalents of the Darcian-scale hydraulic conductivity, in that they have the same dimensions [LT −1 ] and relate a flux at a particular time to a difference in average hydraulic heads at that time between connected but separate bodies of water.They depend on the geohydrological parameters that characterize the subsurface, on the IC and BC (H 0 and H A ), on the forcing parameters R 1 and R 2 , and on time.The Darcian property that the magnitude of the flux does not affect the hydraulic conductivity is maintained in the upscaled conductivity, but the upscaled conductivity no longer is purely defined by properties of the porous medium and the fluid: external forcings also affect it.In this sense, saturated flow at the field scale is fundamentally non-Darcian, even for uniform media and uncomplicated flow patterns.The dependence on R 1 and R 2 implies that the upscaled conductivities change abruptly at t 1 , when the recharge rate changes instantaneously.The time dependence makes them also change gradually.This dependency upon the forcings at this scale arises directly from the expressions for the flux and the average hydraulic head, and exists despite the absence of heterogeneity.The repercussions of this fundamental non-Darcianity at this scale for large-scale groundwater flow modelling with model cells much larger than the Darcian scale are not entirely clear, although it seems reasonable to expect that the assumption of Darcianity at super-Darcian scales becomes increasingly compromised as the curvature of H (x * ,t * ) with x * within the volume of interest increases.All this results in a non-unique Q−H relationship that is more a by-product of the full solution than a tool for flow calculations.The three-variable relationship between the flux across the system boundary, average hydraulic head, and time is unique, but of little practical interest: it is defined for a particular configuration of system parameters and thus changes when the initial and boundary conditions and/or the recharge forcing are changed, even if the geohydrology remains the same.Furthermore, if the problem is fully defined, not only the upscaled hydraulic conductivities can be computed directly, but also the flux as a function of time, through Eqs. ( 38) and ( 47), which do not require H (x,t).Still, Eqs. ( 42)-( 43) and ( 52)-( 53) provide a rigorously derived, insightful, and surprisingly direct relationship between two field-scale subsurface flow characteristics.

Materials and methods
The solutions were coded in Excel worksheets (available upon request from the author), allowing maximum portability and giving considerable flexibility in selecting the values of x, r, and t for output.After some testing, all infinite sums were calculated with 2001 terms.For most values of x, r, and t, this was much more than needed, with the last 100 or even 1000 terms adding less than 10 −3 to the total sum.Only immediately after a change in boundary conditions and x or r equal to 0.99 or 1 did errors up to a few percent remain.In those cases, the error with 1001 terms in the sum was not much larger than that with 2001 terms.Calculation times on a standard laptop ran from too small to determine to about 10 s.
The solutions were used to run a variety of scenarios involving the nine cases discussed above to study the effect of flow geometry (parallel vs. radial), recharge/loss independent of the hydraulic head, and exchange with a lower aquifer.The scale of the modelled problems roughly represents that of an agricultural field with artificial drainage by ditches or tube drains.The labels of the scenarios and the corresponding parameter values are given in Table 1.For the reference cases (Table 1) and the case of 1-day recharge, the accuracy of using truncated series (only the leading term or the first five terms) of the various series appearing in the equations for H (t), Q(0,t), K up , and K ur was evaluated.
Given the somewhat ambiguous character of the storage coefficient (see Appendix A), the sensitivity of flux across the groundwater-surface water interface to variations in µ was examined in some detail.For leaky aquifers, H 0 < H A , while H 2 was set to exceed H A , causing the flux across the groundwater-surface water interface to switch from lateral infiltration to discharge.The effect of the temporal distri- (Table 1) in an infinite strip-aquifer with parallel flow.3 Fig. 2. The evolution in space and time of the hydraulic head for the reference problem (Table 1) in an infinite strip-aquifer with parallel flow.
bution of recharge was studied by providing similar amounts of recharge either uniformly distributed over the first day of the simulation period, or as a pulse at the beginning or the end of that first day.

General
The hydraulic head as a function of space and time for the parallel reference problem (Table 1) is given in Fig. 2. The initial wetting from the ditch during the initial period without rainfall clearly shows, until the system is essentially at equilibrium at t = 0.60 (after 40 days).Despite the logarithmic time scale, it is still clear that the subsequent small recharge flux leads to steady-state flow at t = 2.12, within 41 days after the start of recharge.The plot for the analogous radial problem (Fig. 3) illustrates how the gradients in radial flow can be much smaller.The effect of flow entering/leaving the system in all radial directions leads to much quicker equilibrium than does flow in a singular direction: at t = 0.38 (25 days) is the aquifer at equilibrium with the surface water level, and at t = 1.89 (26 days after the start of the recharge), the flow is steady (compare Figs. 2 and 3).
Under favourable circumstances, the series converge rapidly enough for their first terms to provide accurate approximations (but it is recommended to verify this).In that case, the aquifer behaves like a linear reservoir, with the reservoir coefficient determined by K, D, L (for strip aquifers only), and the aquifer geometry (strip or circular).The full solutions are more flexible than linear-reservoir approximations in that they can handle leaky aquifers, a more flexible set of forcings, and lateral flow towards the surface water as well as lateral inflow from the surface water.This may help explore the deviations from linear-reservoir  1.5 1.5 1.5 behaviour of groundwater reservoirs discussed by Fenicia et al. (2006).Another advantage is that they provide a full map of the hydraulic head in the space-time domain, allowing the results to be compared to, and possibly calibrated on, monitoring well data.

Accuracy of truncated series
If sufficiently accurate results can be obtained with only a few terms of the series, efficiency can be increased and implementing the solutions becomes easier.Since H is expected to be dominated by low-frequency terms, accurate approximations based on the first few terms of the series may well be possible.The flux at the groundwater-surface water interface, on the other hand, is determined by the local gradient in H , which is affected by terms of any frequency, and a (much) higher number of terms may be necessary.This should then also result in loss of accuracy in the estimates of K up and K ur based on truncated series.
Indeed, Fig. 4 (A for the reference case and B for the case with one-day recharge) show that even a single term approximation of H works well except very shortly after an abrupt change of the surface water level.The flows covered by these figures involve flows into and out of the surface water, and are driven by water level changes in the surface water as well as by groundwater recharge by rainfall.Shortly after a change in H A , the fluxes require at least five terms of the series to avoid massive errors (Fig. 5a, early times).The accuracy eventually becomes excellent even for a single term (Table 1) in a circular aquifer with radial flow.3 Fig. 3.The evolution in space and time of the hydraulic head for the reference problem (Table 1) in a circular aquifer with radial flow.) for the reference problem (A; Table 1) and the problem with one-day rainfall (B; Table 1: even rain) in an infinite strip-aquifer (par.) and a circular aquifer (rad.).
(indicating linear-reservoir-type behaviour), but as soon as the flux becomes driven by recharge (after 100 d), the singleterm approximation fails.In Fig. 5b, the flux during the first day is driven by rainfall-generated recharge, and clearly both truncated series underestimate the true discharge flux.During the drying period after the cessation of rainfall, the accuracy rapidly improves.Still, this can give rise to serious mass balance errors over a given period: more water was added by the rainfall than would be eventually discharged if an insufficient number of terms is computed.In cases where five terms are insufficient, often between 1000 and 2000 terms are needed.Still, computational demand (10 s or less on a standard laptop) was not an issue (see Sect. 3). ) for the reference problem (A; Table 1) and the problem with one-day rainfall (B; 4 Table 1: even rain) in an infinite strip-aquifer (par.) and a circular aquifer (rad.).The flux in 5 the infinite strip-aquifer is indicated on the left vertical axis, that in the circular aquifer on the 6 right vertical axis.7 ) for the reference problem (A; Table 1) and the problem with one-day rainfall (B; Table 1: even rain) in an infinite strip-aquifer (par.) and a circular aquifer (rad.).The flux in the infinite strip-aquifer is indicated on the left vertical axis, that in the circular aquifer on the right vertical axis.
The upscaled hydraulic conductivities in Fig. 6 reflect the deviations of the flux for truncated series.Using only the first term results in constant values, consistent with the observation in the Theory section that using only the first terms of the series occurring in the expressions for the relationship between H and Q(0, t) linearizes the relationship.Figure 6 demonstrates the penalty for this simplification: particularly after sudden changes and during periods of recharge-driven fluxes, the single-term approximation fails.For rechargedriven fluxes, even five terms are not enough, as seen by the deviations in Fig. 6a after 100 d and in Fig. 6b during the first day.The results for the leaky aquifers and for the scenario with the smallest storage coefficient (0.001) showed no major deviations from these findings.

Sensitivity to the storage coefficient
The solutions for parallel (Eq.28) and radial flow (Eq.35) show that the storage coefficient appears in exponents that contain the form −t/µ and thus acts as a scaling factor for time: small values of µ make the system respond faster.The smallest values of µ (0.001 and 0.01) reflect conditions for confined aquifers, while the values from 0.1 to 0.4 are more representative for phreatic aquifers.
The anticipated slowing down of the aquifer response is illustrated by the fluxes from the surface water into the aquifer in Fig. 7.The aquifers with large µ not only respond more slowly, but also require more water to fill the larger volume of available storage when water infiltrates laterally from the surface water (until day 100).The fact that the amount of storage in the circular aquifer diminishes farther away from  1) and the problem with one-day rainfall (B; Table 1: even rain) in an 4 infinite strip-aquifer (par.; Kup) and a circular aquifer (rad.; Kur). 5 Fig. 6.The quality of approximation of the upscaled hydraulic conductivity by using 1 or 5 terms of the series solutions compared to the solutions with 2001 terms (ref. ) for the reference problem (A; Table 1) and the problem with one-day rainfall (B; Table 1: even rain) in an infinite strip-aquifer (par.;K up ) and a circular aquifer (rad.;K ur ). the surface water is reflected in its faster response compared to the infinite-strip aquifer.
After recharge starts at day 100, the hydraulic head rises more swiftly in aquifers with small µ, and consequently, the flux to the surface water increases more rapidly (Fig. 7).Here too, the circular aquifers reach near steady-state flow more rapidly than the strip aquifers.with parallel flow that receives water from a deeper aquifer (Table 1: leaky).3 Fig. 8.The evolution in space and time of the hydraulic head in an infinite strip-aquifer with parallel flow that receives water from a deeper aquifer (Table 1: leaky).
The response to variations in µ is marked.While the correct value of µ may be difficult to establish a priori, the sensitivity of key model output to its value makes it a suitable calibration parameter.

Leaky aquifers
The values of parameters a and b for the case of a leaky aquifer are consistent with a head H 2 in the lower aquifer of 4.0 m and a resistance of the aquitard between the aquifers of 100 d.With the imposed BC of 1.5 m and the initial head of 1.0 m, this leads to an appreciable flux into the top aquifer that needs to be discharged into the surface water.The hydraulic head in the top aquifer rapidly rises, and the H -profile slopes down towards the surface water to facilitate the lateral discharge flux (Fig. 8).At t = 0.38 (25 d) for both parallel (Fig. 8) and radial flow (Fig. 9), the hydraulic head has become nearly steady and is entirely determined by the heads in the surface water and in the lower aquifer.The additional head-independent recharge by infiltrating rain of 5 mm d −1 after t = 1.5 (100 d) is small compared to the recharge from below and causes only a minor increase in H .The effect of the dimensionality of the flow manifests itself predominantly in the larger gradients required by parallel flow to drive the lateral flow, resulting in larger deviations from H A overall.This is confirmed by Fig. 10: the response of H to the forcings is more pronounced for parallel flow.During the early stage, where water moves in from the surface water (negative fluxes, owing to the fact that H 0 < H A ), its value is lower for the infinite strip aquifer compared to the radial aquifer.Later, when the influx from the deeper aquifer needs to be laterally transported to the surface water, the infinite strip aquifer has the highest values.The initially negative fluxes gradually trend to zero and become positive within 2.5 d for both 1 Figure 9.The evolution in space and time of the hydraulic head in a circular aquifer with 2 radial flow that receives water from a deeper aquifer (Table 1: leaky).3 Fig. 9.The evolution in space and time of the hydraulic head in a circular aquifer with radial flow that receives water from a deeper aquifer (Table 1: leaky).parallel flow and radial flow (Fig. 11).In the mean time, H gradually increases from being smaller than H A to exceeding it (shortly after 2.5 d for parallel flow and shortly before 2.5 d for radial flow).Thus, there is a brief period during which the direction of the local flux at the groundwater-surface water interface is inconsistent with the magnitude of the field-scale H with respect to H A : the flux is in the direction of increasing average H . Consequently, the field-scale hydraulic conductivity becomes negative during this period (Fig. 10, only visible for parallel flow).
The fluxes infiltrating into the aquifers early in the simulations are comparable (Fig. 11).The discharge fluxes generated by the influx from the lower aquifer and, after 100 d, by recharge from above, are markedly different for the infinite strip and the radial aquifer, reflecting their geometries.

Recharge distribution in time
Three recharge regimes were tested, all for rain showers delivering 0.02 m of water to the aquifer: two involved a pulse application at t = 0 d or t = 1 d (dimensionless time 0.015), in the third regime water entered the aquifer evenly distributed over a one-day period starting at t = 0 (Table 1, problem labels "parallel/radial pulse 0", ". . .pulse 1", and ". . .even", respectively).The pulsed applications raised the hydraulic head by 0.10 m (0.02/µ), followed by a decay back to H A (1.5 m).During evenly distributed recharge, the peak hydraulic head was obviously reached after 1 d: 1.58 m for parallel flow and 1.56 m for radial flow.Within three days, the difference in average hydraulic head was about 1 cm for both parallel and radial flow (Fig. 12).The average head in the circular aquifer dropped at more than twice the rate of the strip aquifer.
The fluxes towards the surface water generated by the recharge show comparable behaviour (Fig. 13).The fluxes for the strip aquifer and the circular aquifer differ by a factor 62.8 (2πL) owing to the necessity of expressing the flux from the strip aquifer per unit length.For the pulse applications, the infinite head gradient at the groundwater-surface water interface caused by the spiked increase in H makes the initial flux go to infinity.The peak discharge for evenly distributed recharge is 0.062 m 3 d −1 per meter for the strip aquifer and 3.4 m 3 d −1 for the circular aquifer.After 3 days, the difference between the largest and the smallest flux has already dropped below 9 % of its peak for the strip aquifer and below 12 % for the circular aquifer.From then on, the flux decays exponentially, with the difference between the three rainfall regimes decaying exponentially as well (the distance on the log scale remains about constant).After 20 days, the 66 1 Figure 12.The average hydraulic head during and after 0.02 m recharge delivered as a pulse 2 at the start (0) or the end of day 1 (1), or evenly distributed over the day (0-1), for an infinite 3 strip-aquifer (p) and a circular aquifer (r).The input parameters are given in Table 1 (pulse 0, 4 pulse 1, even rain).5 Fig. 12.The average hydraulic head during and after 0.02 m recharge delivered as a pulse at the start (0) or the end of day 1 (1), or evenly distributed over the day (0-1), for an infinite strip-aquifer (p) and a circular aquifer (r).The input parameters are given in Table 1 (pulse 0, pulse 1, even rain).flux from the strip aquifer is less than 1.5 % of the peak for any of the three regimes.The circular aquifer loses its water much faster: after ten days, all fluxes are around 1 % (1.1 % at most) of the peak.The conclusion is warranted that for phreatic aquifers, the effect of the temporal distribution of recharge becomes negligible within a few days.For most purposes, daily sums of net infiltration into the saturated zone will likely suffice on input.Since the PDE was linearized, the solution for more complicated rainfall regimes can be obtained by superimposing the solutions for a sequence of daily inputs.

Summary and conclusions
Solutions for linearized parallel and radial flow in aquifers were developed that have sufficient generality to be directly applicable to a wide range of forcings that cover most conditions occurring in nature.The solutions are valid for discharge as well as lateral infiltration.Expressions for the flux between groundwater and surface water, and for the average hydraulic head, were developed.While the test cases were geared towards artificially drained fields, the solutions can be readily applied to hydrological systems with a much sparser discharge network.The linearity of the PDEs allows the possibility to represent the forcing in a given time period as a sequence of time segments (possibly daily segments), the solutions of which can be superimposed to acquire the solution for longer time periods with varying boundary conditions and head-independent and head-dependent recharge.To do so, H 0 for a subsequent solution can be made equal to H resulting from the solutions for previous time periods.This introduces a small error that is likely to damp out rapidly.
The solutions take the form of infinite series, but truncated series of five terms or even a single term can still give accurate results, particularly at times without exchange of water with another aquifer or the unsaturated zone above (if the aquifer is phreatic).When there is such an exchange (headindependent and/or head-dependent recharge), or shortly after a change in the boundary conditions, truncated series will lead to significant errors.
The expressions for the average hydraulic head, the flux across the system boundary (to/from the surface water), and the field-scale hydraulic conductivity contain infinite series.The Darcian nature of the upscaled hydraulic conductivity for steady-state flows can be preserved for transient flows by truncating the series appearing in the expressions for the upscaled conductivity after the leading term.This causes the aquifer to behave like a linear reservoir, with the reservoir coefficient reflecting porous media properties, and the dimensions and geometry (strip/circular) of the aquifer.Using only the first term can lead to significant errors.It is not yet clear how the fundamental non-Darcianity at super-Darcian scales that the solutions prove to exist even in homogeneous media affects large-scale groundwater flow modelling.
The results presented here demonstrate that even for aquifers with small lateral flow distances, a temporal resolution of recharge of one day will often be sufficient.This suggests that soil water models should primarily focus on the accurate partitioning of the precipitation flux into evapotranspiration, direct delivery to the surface water via flow routes that bypass the groundwater (e.g.surface runoff and flow through crack networks and other macropores), and groundwater recharge.Adequate discharge estimates require the model for groundwater-generated discharge described here to be supplemented with a surface runoff and bypass flow model to capture the rapid discharge generation.

Appendix A
The magnitude of the storage coefficient For confined aquifers, the storage coefficient µ reflects the compressibility and consolidation effects on the matrix (the increase in the amount of water stored across D with an increase of the hydraulic head) and generally lies between 5 × 10 −5 and 5 × 10 −3 to 0.01 (Bouwer, 1978, p. 31;Kruseman and de Ridder, 1990, p. 23).For phreatic aquifers, it reflects the storage change across the vertical cross section with a change in the hydraulic head: where x 3 (L) is the vertical coordinate (set to zero at the top of the impermeable layer), θ is the volumetric water content, and x a (L) is the height of the soil surface above the impermeable layer.If instantaneous hydraulic equilibrium is assumed in the vertical (consistent with the Dupuit assumptions), the following equality holds if the phreatic level x p [L] is sufficiently deep for the water content at the soil surface to be equal to θ r : where h is the matric potential [L].This equality can be derived from a plot of the water content against x 3 by first integrating θ over the range of x 3 and subsequently integrating the depth range for each water content over the range of θ.
The integral on the RHS is a soil-specific constant.If H is changed, the Dupuit assumptions stipulate that this results in a similar change in x p .With Leibniz' rule (Abramowitz and Stegun, 1965, 3.3.7),the resulting storage change in Eq. (A1) can be found: In the other extreme, the capillary fringe reaches to the soil surface.In that case, the storage change is similar to that of a confined aquifer.Thus, for phreatic aquifers, µ ranges from nearly zero to (θ sθ r ).Kruseman and de Ridder (1990, p. 23) give a range between 0.01 and 0.30.For conductive aquifers with significant horizontal flow (low clay content), often 0.25 < (θ sθ r ) < 0.4 (compare Table 1.3 of Kruseman and de Ridder, 1990, p. 24).Therefore, µ ≈ 0.2 may be a reasonable approximation.

Appendix B
The solution to the parallel flow problem The PDE for Case 9 is of the form (compare Eq. 5): where the definitions of A and B follow from Eq. ( 5).The IC and BC are: This is a parabolic, non-homogeneous 2nd order PDE with non-homogeneous BCs.The following substitution removes the term AH (see also Farlow, 1993, pp. 58-61): The system becomes: The non-homogeneous term in the PDE is now independent of H .The BC at x = 1 is non-homogeneous.We seek a substitution that makes both BCs homogeneous (compare Farlow, 1993, p. 43-47).Thus: To ensure homogeneous BCs for V (x,t), β(t) needs to be defined as: From Eqs. (B11) and (B13) follows the definition of ζ (t): The substitution in Eq. (B10) thus becomes: The system now becomes: (Note that the term with dB/dt in Eq. ( B6) is cancelled by a similar term that arises when Eq. ( B15) is differentiated with respect to time.) This system can be solved by the method of eigenfunction expansions (Farlow, 1993, p. 64-70).We start from the solution of the associated homogeneous problem by separation of variables: where X and T are as yet unknown functions.From this follows (Farlow, 1993, p. 33-41): The BCs for V give for X: The general solution for X is: with C and E unknown constants.From Eq. (B22) follows that C = 0. From Eq. (B23) then follows that cos(λx) = 0. Thus, we have: Equation (B25) gives a valid solution for arbitrary values of E n .These are therefore set to 1.The non-homogeneous term in Eq. (B16) needs to be expressed as a series of X n : Note that X m and X n are orthogonal on the integration interval for m = n (see Bruggeman, 1999, p. 701-703), and the integral of their product therefore equals 0. In the sum, only the term for m = n is non-zero: Note that setting the value of E m equal to one above is accommodated by the freedom to determine f m .Evaluating the integrals leads to: With this, the series expansion of the non-homogeneous term is: The solution V (x,t) was assumed of the form: The linearity of the system of Eqs.(B16)-(B19) and the homogeneity of the BCs ensures that any linear combination of solutions is also a solution: T n needs to be chosen such that the following equality holds for an individual solution V n : Including in Eq. (B22) the SOV solution to V to gives: We replace the second spatial derivative of X n by −λ 2 n X n (Eq.B21).This allows X n to be divided out: This is an ordinary differential equation in T n .To simplify the notation we introduce and to simplify Eq. (B34) to: The solution to Eq. ( B37) is: with τ [T] an integration variable, and F n a constant that needs to be determined from the IC.The integral in Eq. (B38) can be written as: From Eqs. ( 5), ( 27), and (B1) follows: where u(t) is the Heaviside step function.With this, the final integral in Eq. (B39) can be written as: where the multiplication with the Heaviside step function ensures that the last term is zero for 0 < t < t 1 .
Inserting this into Eq.(B39) gives Inserting this into Eq.(B38) and taking into account Eqs.(B35) and (B36) gives: For t = 0, this reduces to To find F n we therefore need to expand the IC of V (x,t) in a series of the eigenfunctions X n (x).
And finally, with Eq. (B5) we obtain the solution for H (x,t): where A is given by (compare Eqs. 5 and B1):

G. H. de Rooij: Transient flow between aquifers and surface water
To solve the above system, we seek a transformation that makes the inhomogeneous terms independent of H .The transformation of Eq. (B5) creates problems in the term with 1/r, however.Successive application of the Laplace transform and the finite Hankel transform is another strategy (Bruggeman, 1999, p. 744-748).By doing so, the initial condition is transformed differently from the boundary condition though, and there is no term with (H A − H 0 ), which makes the solution somewhat difficult to interpret.We therefore first introduce a substitution to make the BC homogeneous: The PDE then becomes: with IC and BC: Equations ( 26) and ( 27) for a, b, and R(t) remain unchanged, while the BC at r = 0 (Eq.C2) holds for both H and U .
Eliminating the time coordinate by the Laplace transform (Bruggeman, 1999, p. 652-653) gives: where s is the Laplace variable, the subscript L indicates a transformed variable, and B L (s) follows from the definition of the Laplace transform: The transforms were taken from Abramowitz and Stegun (1965, p. 1020-1029).The Laplace transform removes the time coordinate and therefore the need for an initial condition.Next, the system is simplified further by the finite Hankel transform (Bruggeman, 1999, p. 706-707): with α n the roots of The usefulness of the Hankel transform arises from its transformation of the derivatives with respect to r * or r that appear in the flow equations for radial flows (Bruggeman, 1999, p. 707).For r, the derivatives are transformed as: In the above equations, f denotes an arbitrary function of r, J i (y) is the Bessel function of the first kind and order i, and the subscript H denotes a Hankel-transformed variable or function.The boundary condition at r = 1 (r * = L) as well as conditions imposed at r = 0 are incorporated in the transform defined in Eq. (C15).Transforming Eq. (C9) gives: With the relationship for derivatives of Bessel functions of integer order (Abramowitz and Stegun, 1965, 9.1.30;Carslaw and Jaeger, 1959, p. 198, Eq. 5), the integral in Eq. (C16) can be expressed as: 1 0 rJ 0 (α n r)dr = 1 α n rJ 1 (α n r) Solving Eq. (C16) for U LH and using Eq.(C18) gives α n s + With this, the expression for U LH (s,n) becomes This is the solution in the s-n domain.We now seek the inverse transform from the Laplace domain to the time domain.

1Figure 1 .Fig. 1 .
Figure 1.Sketch of the subsurface structures and its model schematization.The variables 2 defined in the main text.3

Figure 2 .
Figure 2. The evolution in space and time of the hydraulic head for the reference problem 2

Figure 3 .
Figure 3.The evolution in space and time of the hydraulic head for the reference problem 2

Figure 4 . 6 Fig. 4 .
Figure 4.The quality of approximation of the average hydraulic head H by using 1 or 5 terms 2 of the series solutions compared to the solutions with 2001 terms (ref.) for the reference 3 problem (A; Table1) and the problem with one-day rainfall (B; Table1: even rain) in an 4 infinite strip-aquifer (par.) and a circular aquifer (rad.) 5 6

Figure 5 .
Figure 5.The quality of approximation of the flux across the groundwater -surface water 2 interface by using 1 or 5 terms of the series solutions compared to the solutions with 2001 3 terms (ref.) for the reference problem (A;Table 1) and the problem with one-day rainfall (B; 4

Fig. 5 .
Fig. 5.The quality of approximation of the flux across the groundwater -surface water interface by using 1 or 5 terms of the series solutions compared to the solutions with 2001 terms (ref.) for the reference problem (A; Table1) and the problem with one-day rainfall (B; Table1: even rain) in an infinite strip-aquifer (par.) and a circular aquifer (rad.).The flux in the infinite strip-aquifer is indicated on the left vertical axis, that in the circular aquifer on the right vertical axis.

Figure 6 .
Figure 6.The quality of approximation of the upscaled hydraulic conductivity by using 1 or 5 2 terms of the series solutions compared to the solutions with 2001 terms (ref.) for the reference 3 problem (A;Table 1) and the problem with one-day rainfall (B; Table 1: even rain) in an 4

Figure 7 .Fig. 7 .
Figure 7. Sensitivity of the flux across the groundwater -surface water interface to the value 2 of the storage coefficient μ in the infinite strip aquifer (A) and the circular aquifer (B).The 3 parameter values for all cases are given in Table 1 (parallel/radial mu ..). 4 Fig. 7. Sensitivity of the flux across the groundwater -surface water interface to the value of the storage coefficient µ in the infinite strip aquifer (A) and the circular aquifer (B).The parameter values for all cases are given in Table 1 (parallel/radial µ ..).

Figure 8 .
Figure 8.The evolution in space and time of the hydraulic head in an infinite strip-aquifer 2

Figure 10 .Fig. 10 .
Figure 10.The average hydraulic head and the upscaled hydraulic conductivities (Kup, Kur) for 2 parallel (p) and radial flow (r) in leaky aquifers (see Figs. 8 and 9) with head-dependent 3 recharge from a deeper aquifer.4 Fig. 10.The average hydraulic head and the upscaled hydraulic conductivities (K up , K ur ) for parallel (p) and radial flow (r) in leaky aquifers (see Figs. 8 and 9) with head-dependent recharge from a deeper aquifer.

Fig. 11 .
Figure 11.As Fig. 10, for the fluxes across the groundwater -surface water interface.Th 2 flux in the infinite strip-aquifer (par.) is indicated on the left vertical axis, that in the circula 3 aquifer (rad.) on the right vertical axis.4 Fig. 11.As Fig. 10, for the fluxes across the groundwater -surface water interface.The flux in the infinite strip-aquifer (par.) is indicated on the left vertical axis, that in the circular aquifer (rad.) on the right vertical axis.

5 Fig. 13 .
Figure 13.As Fig. 12, but for the fluxes across the groundwater -surface water interface.The 2 flux in the infinite strip-aquifer (p) is indicated on the left vertical axis, that in the circular 3 aquifer (r) on the right vertical axis.4 5 Fig. 13.As Fig. 12, but for the fluxes across the groundwatersurface water interface.The flux in the infinite strip-aquifer (p) is indicated on the left vertical axis, that in the circular aquifer (r) on the right vertical axis.

Table 1 .
Parameter values for various illustrative problems.

Table D1 .
List of symbols.