A critical assessment of simple recharge models : application to the UK Chalk

Quantification of the timing and magnitude of point-scale groundwater recharge is challenging, but possi- ble at specific sites given sufficient high spatial and temporal resolution field observations, and a suitable physically based model. Such models are generally too computationally in- tensive and have too many unknown parameters to be practi- cally applicable within distributed, larger-scale hydrological or groundwater models. This motivates the need for simpler recharge models, which are widely used within groundwa- ter models. However, it is important that these models are able to capture adequately the unsaturated zone flow pro- cesses. We perform an inter-comparison of recharge sim- ulated by a detailed physically based model and a simple recharge model, with both models applied to a field site in the fractured porous Chalk in the UK. Flow processes are simu- lated convincingly using a dual permeability, equivalent con- tinuum, vertically heterogeneous, Richards' equation model, applied to a 2-D hillslope transect. A simple conventional recharge model was then calibrated to reproduce the wa- ter table response simulated by the physically based model. The performance in reproducing the water table was surpris- ingly good, given the known discrepancies between the ac- tual processes and the model representation. However, com- parisons of recharge fluxes simulated by each model high- lighted problems with the process representations in the sim- ple model. Specifically, bypass flow events during the sum- mer were compensating for recharge that should have come from slow, continual drainage of the unsaturated zone. Such a model may still be useful for assessment of groundwater resources on a monthly basis, under non-extreme climatic conditions. However, under extreme wet or dry conditions, or under a changed climate the predictive capacity of such models is likely to be inadequate.


Introduction
Quantifying the timing and magnitude of groundwater recharge, that is the flux of water across the water table into an unconfined aquifer, remains a significant challenge because there is no direct, non-destructive method for measuring this flux in the field.Indirect methods include using chemical measurements (i.e.tracers) and/or physical Introduction

Conclusions References
Tables Figures

Back Close
Full measurements (e.g.water content, water table response) (Scanlon et al., 2002).With all of these methods, the estimation of recharge requires a conceptual model (that is, a set of assumptions about how the system functions) formulated as a mathematical model (from which quantitative recharge estimates are calculated).One of the problems with using tracers, especially in dual permeability systems, is that the same chemical profile could be consistent with multiple conceptualisations (a classic example is the reinterpretation of the results of Smith et al. (1970), by Foster (1975), discussed in Mathias et al. (2005)).Methods based on physical measurements include: soil water balances (Rushton, 2005;Ragab et al., 1997;Finch, 1998); solutions to some form of 1-D Richards' equation involving soil moisture and/or matric potential observations (van den Daele et al., 2007;Brouy ère, 2006;Habets et al., 2010;Ireson et al., 2009b;Ireson and Butler, 2011); and water table fluctuations (Cuthbert, 2010;Ireson et al., 2012).All of these methods have limitations, which may or may not be prohibitive for a given field site.Soil moisture calculations do not rigorously account for the delay and attenuation of the soil drainage by the unsaturated zone, which can be substantial in aquifers with large unsaturated zones (Rushton, 2005).1-D Richards' equation approaches cannot rigorously represent unsaturated-saturated zone interactions (Ireson et al., 2009b;Ireson and Butler, 2011), which may be important when water table fluctuations cause significant changes in the unsaturated thickness.As well as the specific yield, water table fluctuation methods require estimates of the net groundwater drainage, typically made using a simplified conceptualisation of the saturated zone (e.g.Fig. 1 in Cuthbert, 2010).Therefore, when addressing the problem of quantifying recharge for any given location, careful consideration of the flow processes in the unsaturated and saturated zones needs to be taken.The most appropriate method will depend on this conceptual understanding of the system (Scanlon et al., 2002).
There is a long history of research into the physical processes that control recharge to the unconfined Chalk aquifers of north west Europe, particularly in the UK (see reviews in Mathias et al., 2006;Ireson et al., 2009b;van den Daele et al., 2007;Butler et al., 2012), France (Habets et al., 2010) and Belguim (Brouy ère et al., 2004;Introduction Conclusions References Tables Figures

Back Close
Full  , 2006).The work in these citations and others (Ireson and Butler, 2011;Ireson et al., 2012) has focused on developing physically based models, and using these models to reproduce field observations.However, these have not yet been implemented into operational distributed groundwater models.The primary reason for this is that these models are too complex, requiring a detailed and expensive numerical solution to a system of non-linear, partial differential equations.This results in a large number of parameters and a time consuming solution.Other workers have focused on the more pragmatic task of developing a conceptual recharge model for the Chalk, with a less strict focus on the detailed unsaturated zone processes, and more of an emphasis on the larger groundwater system response (Rushton and Ward, 1979;Howard and Lloyd, 1979;Morel, 1980;Jackson and Rushton, 1987;Finch, 1998;Ragab et al., 1997;Finch, 2001;Bradford et al., 2002;Rushton, 2005).This paper seeks to explore what insights a detailed model can provide about the performance of simpler, conventional recharge models, which have not been so rigorously tested against field observations from the unsaturated zone.Rushton (2005) provides a good, authoritative summary of conventional approaches for modelling recharge to Chalk aquifers.These typically involve three components.(1) A near surface water balance calculation (such as the Penman (1950) and Grindley (1967) approach, see also Sect.8 of Allen et al. (1998)), which determines soil moisture deficit (SMD), actual evaporation (AE) and soil drainage.(2) A fixed proportion of rainfall and/or rainfall over some threshold (explored in some detail by Rushton and Ward, 1979) which becomes bypass flow that passes directly to the aquifer through the fractures (this dates back to the widely cited study by Smith et al. (1970) and was first applied by Mander (1977Mander ( , 1978))).( 3) Where the unsaturated zone is deep (i.e.≥ 10 m) delayed recharge can occur.Rushton (2005) comments that for "chalk and limestone aquifers the response in observation boreholes is usually within one month of the occurrence of recharge", with exceptions associated with very deep unsaturated zones containing marl bands (see for example Cross et al., 1995).Findings from previous studies (Ireson et al., 2009b;Ireson and Butler, 2011) would imply that even for a relatively uniform Chalk profile, each of these components Introduction

Conclusions References
Tables Figures

Back Close
Full is more complicated.In this paper we will try to address these complications explicitly by comparing detailed, physically based models calibrated against field observations with simpler, conventional recharge models.

Method
The physically based model for 1-D vertical flow in the Chalk unsaturated zone presented by Ireson et al. (2009b) is extended to represent 2-D vertical and horizontal flow in an idealized hillslope transect.Richards' equation governs flow in the saturated and unsaturated zones, and is solved numerically.The properties of the Chalk fractures and matrix are represented using an equivalent continuum approach (Peters and Klavetter, 1988;Doughty, 1999;Ireson et al., 2009b).The original 1-D model was driven using observations of precipitation and actual evaporation, calibrated against near surface observations of water content and matric potential, and used to simulate deep recharge fluxes.As lateral flows are not explicitly represented in the 1-D model a fixed water table boundary condition was used at the lower boundary (40 m and 75 m below ground level for the two field sites).In this paper, by simulating horizontal flow and, in particular, including the actively flowing portion of the saturated zone within the model, we are able to simulate explicitly the observable water table response, rather than the unobservable recharge fluxes alone.This, in principle, provides an opportunity to test in a more rigorous manner the recharge model against observations.However, this is not without complications.Firstly, since the resulting 2-D dual-permeability numerical model is very computationally demanding, it is not feasible to search the parameter space effectively.Secondly, in a real catchment, ideal, fixed transects do not exist -the interfluve (i.e.no flow boundary) moves as the groundwater catchment expands and contracts (a particular issue in the Chalk as shown by Parker, 2011), flows converge onto and diverge away from apparent flowlines (Troch et al., 2003) and flow directions change.One solution to this latter problem would be to move to a 3-D model.However, this is significantly more challenging since, in addition to the added computational expense, in three-dimensional natural systems it is difficult to identify Introduction

Conclusions References
Tables Figures

Back Close
Full static domain boundaries for which states (i.e.type 1 boundaries) or fluxes (i.e.type 2 boundaries) can be quantified (a problem that is sometimes addressed in groundwater models by simulating an area that extends beyond the domain of direct interest, see Parker, 2011, Jones et al., 2012).
In this paper we attempt to capture the essential physical behaviour of the water table response to recharge in a fractured porous Chalk catchment.The approach is to first establish a plausible 2-D hillslope transect model (the benchmark model), which is broadly consistent with field observations (Sect.2), and second, to attempt to emulate the behaviour of this model with a simpler, parsimonious recharge model coupled with a groundwater model (Sect.3).This approach allows us to test many of the assumptions within groundwater models with no explicit representation of the unsaturated zone.It also allows us to explore the identifiability of parameters of conventional recharge models, and the associated problem of equifinality (Ebel, 2006) of calibrated groundwater models.
2 Physically based hillslope model

Field site and data
This study focuses on a transect through the unconfined Chalk in the Pang catchment, a tributary of the River Thames in Berkshire, England (Fig. 1).The site is attractive because of the data availability, and the relatively well understood flow paths in the saturated zone (Hughes et al., 2011;Parker, 2011).The Pang and adjacent Lambourn catchments have been extensively monitored over a number of decades.Four monitoring wells are useful for our study, Knollend Down, Malthouse, Hodcott and Compton (Table 1).An additional observation well, East Ilsley, also lies on this transect, and is included for reference purposes since recent studies have focused on this site (Ireson et al., 2009a(Ireson et al., , 2012;;Gallagher et al., 2012).However, no data were available for East Ilsley before 2003.Introduction

Conclusions References
Tables Figures

Back Close
Full The transect of 5 monitoring wells runs along a (normally) dry valley in the Pang catchment.This is assumed to approximate a groundwater flow line (Hughes et al., 2011), with the lowermost piezometer located at Compton, which can be considered the start of the continuously flowing section of the river Pang -i.e. the lowermost part of the intermittent portion of the river under normal conditions.A digital elevation model indicates that the dry valley weaves downward from Knollend Down to Compton.The stratigraphy across the transect was extracted from the British Geological Survey's 3-D geological model, and is consistent with borehole logs.The Chalk subdivisions are based on broad hydrogeological properties and show as the higher permeability middle Chalk and the lower permeability lower Chalk.

Richards' equation model
The hillslope model used in this study is based on extending the model for 1-D vertical flow in the Chalk unsaturated zone presented by Ireson et al. (2009b).Flow processes were simulated using Richards' equation, solved using a two-dimensional finite volume model applied on an unstructured mesh (Narasimhan and Witherspoon, 1976;Bear and Cheng, 2008, Sect. 8).The governing equation is where S is the volume of water (m 3 ), V is the cell volume (m 3 ), A is the area of each face of the cell (m 2 ), θ is the volumetric water content within the cell (-), t is time (d), q is the Darcy flux (m d −1 ) and s is a sink term (d −1 ) representing uptake of water by transpiration (described below).The fluxes are given by Darcy's law, i.e.
Full where ψ is the pressure head (m), K (ψ) is the hydraulic conductivity (m d −1 ), z is elevation above datum (m), and L is the length between two nodes (m).The model domain is discretised into nodes and cells using a Voronoi, or perpendicular bisection, grid (Palagi and Aziz, 1994), whereby cell faces are located midway between adjacent nodes and are perpendicular to the line that bisects the two nodes.Model states (ψ and θ) are specified at nodes, and fluxes are approximated across faces.Using finite difference approximations, the flux between cells m and n, q m,n , is given by The hydraulic conductivity is estimated at the boundaries using the arithmetic mean (Parissopoulos and Wheater, 1988).The integral of all fluxes into cell m is approximated by where f m is the number of faces associated with cell m and A m,n is the cross-sectional area of the face shared by cells m and n.
The sink term on the right hand side of Eq. ( 1) is approximated as In order to solve Richards' equation for both saturated and unsaturated conditions, it is convenient to solve for pressure head rather than water content as the dependent variable.Therefore, the left hand side of Eq. ( 1) can be rewritten, and approximated for Introduction

Conclusions References
Tables Figures

Back Close
Full cell m by where S s is the specific storage coefficient (m −1 ), S e (ψ) is effective saturation (-), and ). Note, this is based on the pressure head form of Richards' equation which Tocci et al. (1997) demonstrated to be mass conservative using the solution procedure described below.Substituting Eqs. ( 4)-( 6) into Eq.( 1) and rearranging we obtain a system of ordinary differential equations given by This system of equations is solved for ψ(t, x, z) numerically in MATLAB with the stiff ordinary differential equation solver ode15 s (Shampine and Reichelt, 1997).This employs an adaptive time grid to minimise numerical errors, and boundary conditions are applied on a specified time step.This is similar to the solution procedure applied by Tocci et al. (1997), and has been used previously in other 1-D problems (Ireson et al., 2009b;Ireson and Butler, 2011).

Hydraulic properties
The hydraulic properties of the Chalk are represented using the same Equivalent Continuum approach applied by Ireson et al. (2009b).The assumption is that local exchanges of water between fractures and matrix are instantaneous, such that they are in pressure equilibrium.Ireson and Butler (2011) and Ireson et al. (2012) showed that Introduction

Conclusions References
Tables Figures

Back Close
Full this assumption is suitable for representing the recharge processes under normal and dry conditions.High intensity rainfall events can lead to preferential flow in the fractures.Ireson et al. (2009a) found 11 such events occurred in a two year record at East Ilsley.Of these, 9 could be considered negligible in terms of the volume of water preferentially recharged, and one was associated with a very rare, extreme rainfall event in July 2007 (Marsh and Hannaford, 2007).They estimated that even for this rare event, preferential recharge accounted for no more than 15 % of the annual total recharge in 2007.Preferential recharge is a concern in terms of contaminant transport, and in terms of possible future flood risk.It is probably unimportant from a water resource/drought perspective, since such events are rare and the total volumes likely to be recharged are a small proportion of the annual total.For the purposes of this paper, we attempt to simulate the period 1970-2000 and ignore preferential recharge.
The model contains 17 parameters and we adopt values from Ireson et al. (2009b, Table 2) obtained from the field site at Warren Farm in the adjacent Lambourn catchment.The only parameter that was modified was the fracture saturated hydraulic con- ).This was set to 4000 m d −1 , such that with a fracture porosity of 1 % the bulk saturated hydraulic conductivity of the rock would be 40 m d −1 .This magnitude is consistent with the saturated hydraulic conductivity used in a previous groundwater model developed for the area (discussed in Parker, 2011;Jackson, 2011;Jackson et al., 2011;Jones et al., 2012), which for the middle Chalk layers varies between 4 and 76 m d −1 .

Boundary conditions
The model domain is enclosed by the four boundaries shown in Fig. 1d.An ideal hillslope for simulating flow processes would be located on a transect through a homogeneous material, along a flow line where flow is perfectly parallel with the transect (i.e.does not converge or diverge), and bounded by a fixed water table divide at the interfluve, and a fixed stage river in the valley (i.e. the classic conceptualisation of Introduction

Conclusions References
Tables Figures

Back Close
Full T óth, 1963).In real catchments, such a configuration is unlikely to be found.Here we adopt idealised boundary conditions, with the hope that, at the point of interest close to the centre of the transect, the flow processes in the saturated zone are a reasonable approximation to reality.In the Pang catchment, the groundwater divide moves, as is evident from observed water table surfaces (Hughes et al., 2011;Parker, 2011), and this is likely to a be a characteristic that is common to Chalk catchments.Therefore, it is not possible to specify a no-flow boundary at some fixed interfluve location.We therefore consider a uniform inflow at the left hand boundary (Knollend Down) which is downslope from the actual (moving) interfluve.Increasing this inflow results in a larger gradient in the water table along the transect.The value of this inflow was established by trial and error, to give an improved fit to the observed water table elevations in the boreholes along the transect.This was the only refinement of the physically based model that was performed in order to improve the fit with the observations.The value of this inflow was 6 m 2 d −1 (note, this is a flow per unit width of the aquifer), which is equivalent to the contribution from an upstream semi-infinite aquifer of length 6 km, recharged uniformly at 1 mm d −1 .This magnitude therefore is not unrealistic, although in reality we would expect variations seasonally and between years as recharge rates vary and the catchment expands and contracts.This boundary condition is implemented as a uniform lateral flux beneath the water table on the left hand boundary using the discretised form of Eq. ( 8) given in Eq. ( 9).
where (m) and the index i (-) refers to all submerged cells on the left hand boundary, from the base of the aquifer (i = 1) to the water table (i = N).Above the water table, a zero flux boundary is applied on the left hand side.The right hand boundary, at Compton, is located at the point where the intermittently dry valley meets the flowing river Pang.A borehole at Compton provides a continuous timeseries of groundwater level fluctuations, and a type one boundary condition can be used.However, we found that simulated water levels in the upstream observation boreholes were insensitive to the transient water table fluctuations at this boundary, and that using a constant head equal to the mean observed level at Compton, h R = 85 m, gave the same result (upstream from the boundary), and the model converged more easily, reducing computational expense and runtime for the ODE solver.
The lower boundary (i.e. the base of the actively flow portion of the aquifer) is approximately located at the interface between the middle and lower Chalk, where the hydraulic conductivity drops significantly.For simplicity, this is treated as an impermeable boundary and a no flow boundary condition is used.
Finally, for the upper boundary at the ground surface spatially uniform infiltration of rainfall is applied as a specified flux boundary.The high saturated hydraulic conductivity at the surface means that there is no need to account for surface ponding or runoff.Rainfall is given by a daily rainfall time series from a local raingauge (Environment Agency gauge 268812 located approximately 6 km from the transect).The vertical distribution of evapotranspiration is assumed to be entirely transpiration losses from the root zone, and is estimated using the Feddes model (Feddes et al., 1976), driven by monthly MORECs (location 159, Thompson et al., 1981) potential evaporation (PE) linearly interpolated onto a daily time step.
The sink term in cell m due to evapotranspiration through plant roots, s m (d −1 ), is obtained using the root-extraction function

Conclusions References
Tables Figures

Back Close
Full where z m (m) is the depth of the centroid of the cell below the ground surface, ) is the potential evapotranspiration rate, f (ψ m ) (-) is the water stress function, g(z m ) (-) is the normalised root density, A s is the surface area over which the evaporation flux was applied and V m is the volume of the cell.The water stress function, ignoring the effect of anaerobiosis (that is suppression of evaporation under saturated conditions), is given by where ψ d (m) and ψ w (m) are matric potential thresholds corresponding to the point below which water stress commences and the wilting point, respectively.Thresholds values of ψ d = −4 m and ψ w = −150 m were used (Feddes et al., 1976).
The root density, ρ r (m m −3 ), is assumed to be exponentially distributed with depth, defined by where L (m) is a parameter representing the depth above which 63 % of roots are present, taken as 0.2 m (Ireson et al., 2009b).It is assumed that the root systems are uniform along the transect.The normalised root density, g(z m ) (-), expresses the proportion of the total root system present within the particular depth horizon enclosed by cell m.Cells in the mesh that was used were deliberately arranged in columns, and we adopt here the indices i and j to refer to cells within a particular column.Hence, in discrete notation, g(z i ) is given by

Conclusions References
Tables Figures

Back Close
Full Equations 10 and 13 can be combined to give and actual evaporation, AE (m d −1 ) can be found from

Initial conditions
An arbitrary initial condition was used.This was based on the steady-state water table response to a constant recharge rate of 1 mm d −1 , and a uniform pressure in the unsaturated zone above the hydrostatic capillary fringe of −2 m.The initial condition applied on 1 January 1970, and the first 5 yr are considered a warm up period to eliminate the impact of this arbitrary initial condition.

Model performance
The performance of the model in reproducing the observed behaviour of the water table is shown in Fig. 2. Overall, the performance at Malthouse is the best.Here, the level of the annual peak and trough of the water table is normally captured quite well, as is the shape of the hydrographs -with some exceptions.The observed response appears to lag the simulated response consistently, which is a concern.At Hodcott and Knollend Down there are problems during the drought of 1976/1977 which might be explained by local pumping that occurred to augment river flows during this period (Morel, 1980).This does not show up at Malthouse since the observations bottom out, as can be seen in 1976/1977 and in the early 1990s.Overall, given the necessary simplifying assumptions in constructing the domain and boundaries of the transect, the fact that model was essentially uncalibrated, and uncertainty in the driving rainfall and

HESSD Introduction Conclusions References
Tables Figures

Back Close
Full evaporation data, this is a plausible physically based model of recharge and water table response in the Chalk.Figure 3 illustrates the simulated recharge fluxes as compared with rainfall and actual evaporation, for the years 1981 and 1982.High intensity and sustained summer rainfall events typically result in either a negligible or very small recharge response.Conversely, in the winter months when the soil and unsaturated zones are wetter and the unsaturated thickness is smaller, rainfall events are converted into significant recharge responses.This shows that the complex and non-linear recharge processes in the Chalk that have been described elsewhere (Ireson et al., 2009a,b) are captured by the model.
Figure 4 shows hydraulic head profiles in the unsaturated and saturated zones throughout a typical year, with the lower plot showing the seasonal water table response.Flowlines are effectively perpendicular to the head contours plotted.The flow direction changes very sharply between the unsaturated and saturated zones, as the fractures, with a high air entry value (i.e.thin capillary fringe) and high hydraulic conductivity, rapidly empty.In the saturated zone flow is always sub-horizontal, always downslope, and the rate can be inferred from the distance between contours.In the unsaturated zone flow directions are more complex.Despite appearances in Fig. 4 due to the scale distortion, lateral fluxes in the unsaturated zone are negligible (around four orders of magnitude smaller than the vertical fluxes).However, the vertical fluxes switch direction depending on when infiltration or evaporation dominate (winter and summer respectively), with zero flux planes developing in the profile, as described by Wellings (1984) and Ireson et al. (2009b).

Low-dimensional recharge model
In this section we report two approaches which were applied to try to emulate the behaviour of the complex physically based model described above.In both cases the premise is that the hillslope domain can be represented as a series of vertical Introduction

Conclusions References
Tables Figures

Back Close
Full unsaturated columns which do not interact with one another, overlying a water table, where flow is simulated using the Boussinesq equation, as in Eq. ( 16).
where S y is specific yield (-), assumed equal to the fracture porosity (1 %), K s is the bulk rock saturated hydraulic conductivity (m d −1 ) (in this case 40 m d −1 ) and R is the time varying recharge flux (m d −1 ).
The premise was tested by extracting the simulated recharge fluxes, R(t, x) (specifically the interpolated vertically downward fluxes just above the water table) from the Richards' equation model, and using these to drive a Boussinesq equation model, with boundary conditions where L is the length of the model domain (m) (i.e.x = L is the left hand boundary).The initial condition in h was identical to that employed in the Richards' equation model for the initial water table elevation.This 1-D model was solved numerically using the same procedure described above, with MATLAB ode15s.As shown in Fig. 2, the performance of the Boussinesq model was nearly identical to the Richards' equation model.This validates the Dupuit assumptions (p.163, Fetter, 1994) for this particular problem, and justifies the premise of the emulation strategy: in order to emulate the behaviour of the 2-D model all we need to do is reproduce the recharge fluxes.However, this does not mean that the recharge fluxes are necessarily independent of the water table.Introduction

Conclusions References
Tables Figures

Back Close
Full

Approach 1: emulation by linear models
Our first attempt to simulate the recharge fluxes was based on the integral balance method suggested by Duffy (1996), who posed the question "Can low dimensional dynamic models of hillslope-scale and catchment-scale flow processes be formulated such that the essential physical behavior of the natural system is preserved?".Once a detailed numerical model has been established, the states within and boundary fluxes between particular sub-volumes of the model domain can be obtained by integration.Duffy reduced the complexity of a hypothetical hillslope model into a single saturated store and a single unsaturated store, in order to simulate total groundwater storage and discharge.Here, we aim to simulate recharge for a single column.From the Richards' equation model, we can extract the recharge flux and the total unsaturated zone storage as a function of time, for a single column, and together with the effective rainfall time series (rainfall minus actual evaporation), these sum to zero.The resulting unsaturated zone storage-recharge relationships (recall recharge is the flux out of the unsaturated zone) are shown in Fig. 5, with recharge partitioned into fracture flow, matrix flow and total flow.For this particular problem there is no clear relationship between the storage and recharge, even on monthly and annual times scales.We also looked at subdividing the unsaturated zone into 2 or 3 stores, and looked at using the water table (i.e.saturated zone storage), rainfall and evaporation rates as additional covariates.
Using multiple linear regression we were unable to find effective relationships for reproducing the recharge without incorporating autoregressive terms, e.g.
When autoregressive terms were included, reasonable regression relationships could be found, but the resulting model was unstable, with errors that accumulate with time.
Whilst simple recharge relationships could not be found using this method, we were able to find a simplified evaporation algorithm.In the Richards' equation model, po- where AE and PE are actual and potential evaporation (m d −1 ), S is the bulk water content (i.e. the depth averaged volumetric water content) in the soil zone (-) taken as being 0.275 m deep, and β are the regression coefficients.Rearranging this equation and substituting in the coefficients, we obtain This relationship is shown inset in Fig. 6, and is in fact similar to the water stress function (involving soil moisture deficit and root constant) in standard soil moisture models (Penman, 1950;Grindley, 1967).This suggests that for calculating actual evaporation from the Chalk, these standard models are adequate, and that no significant benefit arises from using the more complicated Feddes model.

Approach 2: calibration of a conventional recharge model
The second approach that was applied was to attempt to fit a simple conventional model to the benchmark water table simulated by the Richards' equation model.The conventional model comprises a simple recharge model, coupled to the 1-D Boussinesq equation.This is similar to the recharge model presented by Rushton (2005).The states, fluxes and parameters used in this model are given in Table 2.A schematic diagram of the model is given in Fig. 7.
The model is solved as an initial value problem, with an explicit time stepping method, with index j refering to the current time step.An initial soil moisture deficit, M A,1 is specified, and then Eqs. ( 20)-( 25 n T , where n T is the number of time steps. B j = 0, (P j − TH) ≤ 0 BF(P j − TH), (P j − TH) > 0 (20) PWP−RC , RC < M P ,j < PWP 0, M P ,j ≥ PWP ( 22) This model was coupled with the Boussinesq model for the hillslope transect, such that recharge was uniform across the entire transect.This was coded up in FORTRAN using a simple explicit solution scheme, with a time step sufficiently small to give acceptably small truncation errors, and was thus very computationally efficient, allowing us to search the parameter space effectively using a Monte Carlo method.There are four parameters in the recharge model, and 10 4 realisations were used.Stratified parameter sampling was adopted such that one-third of realisations had no bypass (BF = 0, TH = 0), one-third had bypass with no threshold (BF > 0, TH = 0), and onethird had bypass over some threshold (BF > 0, TH > 0).These combinations result in three separate model structures, and the sampling gives an equal weight to each structure.This allows us to explore all of the model configurations considered by Rushton and Ward (1979), who were only able to consider 9 parametric combinations due to Introduction

Conclusions References
Tables Figures

Back Close
Full .The objective function minimized the root mean squared error of the water table simulated by the simple and benchmark models, at a point in the centre of the hillslope transect, in order to minimize the impact of the lateral boundary conditions on the response.

Conventional recharge model results
The performance of the simple model in reproducing the benchmark water table is shown in the first plot in Fig. 8a.Overall, the performance is surprisingly good.The most noticeable weakness is during the drought in the mid 1970s, and dry periods in the late 1980s and early 1990s.However, even during these periods these errors are small, especially compared with errors between the benchmark model and observations.To explore parameter sensitivity, the 10 optimal model realisations are plotted, and are all equally good fits to the benchmark water table (equifinality).Parameter identifiability is shown in the dotty plots in Fig. 8c and optimal parameter ranges that are quoted are based on the range of parameter values for the 10 optimal realisations.The bypass threshold is insensitive and therefore should not be included in the model, since it requires an additional parameter that results in no improvement.The bypass fraction is sensitive, with an optimum value in the range 7.1-10.5 %.The root constant is also relatively sensitive, with an optimal value in the range 420-570 mm.The permanent wilting point is less sensitive, but larger values are subtly better. is very unlikely to ever drain completely.This therefore implies that a depth of Chalk significantly greater than 3 m is supplying water for summer evaporation, which is consistent with the deep zero flux planes hypothesized by Wellings (1984) and simulated to depths of 6-7 m by Ireson et al. (2009b).Simulated actual evaporation tends to be somewhat higher in the summer in the simple model as compared with the benchmark model.In 1976, some of the 10 optimal simple model realisations appear to dry to soil to wilting point, with evaporation rates dropping to zero for some time.To better compare the recharge signals, Fig. 8b shows the distribution of recharge fluxes simulated by each model.It can be clearly seen that whilst the model does a reasonable job of reproducing higher recharge rates (R > 1 mm d −1 ), the simple model (i) overestimates the low recharge rates; and (ii) predicts no recharge (or strictly R < 0.01 mm d −1 ) for about 40 % of the time.In the benchmark model, there is always some small amount of recharge from slow drainage of the unsaturated zone (discussed by Ireson et al., 2009b).However, during the summers, and especially during the droughts of the mid 1970s and early 1990s, the simple model is unable to simulate slow drainage from the soil and unsaturated zone which, in reality, sustains recharge rates and water levels.
Instead, the simple model compensates for this by using bypass rainfall events, which are visible in the recharge signal during these periods.These events are not physically realistic, but they do enable the model to transmit enough water down to the water table to maintain higher summer water levels.The calibration thus tailors the bypass fraction to optimize this summer response.This is clearly demonstrated in Fig. 9 where the 10 optimal model realisations with and without bypass are compared.Without bypass, the model is unable to simulate the water table troughs, especially during drought conditions, and model performance is worse overall.Therefore, whilst this model performs well, the physical mechanisms are not correct.Note that if there is no soil moisture deficit, all rainfall in excess of potential evaporation is passed directly to the water table, so bypass becomes irrelevant.As a result, in the winter when the soils are wet and SMD is frequently zero, the effect of bypass is minimal.Introduction

Conclusions References
Tables Figures

Back Close
Full It is also interesting to look at how the same simple model would perform if applied on a monthly time step.Groundwater models tend to be run on monthly time steps, so this is a more realistic test of how simple recharge models might be applied to real problems.Monthly simulations were run by (i) calculating the recharge daily recharge rate as before; (ii) accumulating the simulated daily recharge rates to monthly recharge rates; and (iii) running the Boussinesq equation model on a monthly time step.The results are shown in Fig. 10 (note that since the recharge calculation is still daily, soil moisture and actual evaporation are unchanged and hence not replotted).There is no significant decline in the performance of the simulated water table, shown in Fig. 10a.However, the simulated recharge time series in Fig. 10a and cumulative distribution in Fig. 10b appear significantly improved.This is because, when the recharge rates are accumulated onto a monthly time step, all of the individual summer bypass events noted previously are lumped together and, again due to the calibration, approximately make up the slow drainage from the more physically realistic model.Unsurprisingly, the optimal parameter ranges and the parameter identifiability are different in the monthly model configuration, as shown in Fig. 10c.Again the bypass threshold can be eliminated as unnecessary.Again the bypass fraction is the most sensitive parameter, but now takes a somewhat lower value, in the range of 4.8-6.1 %.The root constant has an optimal value in the range 1190-1250 mm, but in fact performance does not significantly decline for larger values.Finally, the permanent wilting point is now completely insensitive.

Conclusions
In this paper we outline a framework that can be applied for the rigorous quantification of the timing and magnitude of groundwater recharge, taking a UK Chalk aquifer as a case study.For general application, the procedure is as follows: (i) propose a conceptual model for flow processes in the unsaturated and saturated zones, on the basis of understanding grounded in field observations; (ii) develop a physically based, small Introduction

Conclusions References
Tables Figures

Back Close
Full scale model that represents all of the processes as accurately as possible; (iii) apply the physically based model to a well designed field experiment, where boundary conditions and flow processes are well understood; (iv) emulate the essential behaviour of the successful physically based model, using simpler models, ideally based on existing, operational models.Each stage depends on the success of the previous stage, and stage three in particular is challenging, and the entire focus of many research efforts.The final stage will inform the structure and parameters used in recharge models that can be coupled with groundwater models.Following this procedure is not a trivial exercise, and will not be appropriate in all cases.For example, in data sparse regions where annual recharge estimates are all that is needed (or achievable), this would certainly not be a sensible approach.However, the advantage of this method is that it allows us to explicitly test many of the assumptions on which our recharge model (conceptual and mathematical) is based.
In this study we have applied a physically based model for the Chalk unsaturated and saturated zone based on Richards' equation with dual permeability, equivalent continuum properties, to a hillslope transect.Under certain necessary simplifying assumptions about the model domain, and with minimal calibration or refinement, the model performance was broadly consistent with transient field observations of water levels in monitoring wells.This model is highly computationally demanding, and consequently wider application of this model is impractical for anything other than research purposes.This motivates the need for a simplified model that can emulate the simulated responses.
The first attempt to emulate the model failed.This used an integral balance approach (Duffy, 1996) to extract integrated states and fluxes from the Richards' equation model, and attempted identify simple constitutive relationships between these using multiple linear regression.The only part of the system for which a simple relationship could be found using linear modeling was between actual evaporation, potential evaporation and soil moisture.This resulted in a simple power law relationship that is in fact very Introduction

Conclusions References
Tables Figures

Back Close
Full similar to the water stress function found in many simple soil models (e.g.Penman, 1950;Grindley, 1967;Rushton, 2005).
The second attempt at emulation was based on fitting a simple conventional recharge model (after Rushton, 2005) to the water table simulated by the benchmark, physically based model.The simple model performed surprisingly well.We believe this demonstrates the power of the use of storage-threshold type relationships in hydrological models, which are anyway well established (in particular in the models that involve the soil moisture deficit concept), and produce a response that is completely distinct from linear models.However, on close examination of the recharge signal simulated by the simple model, it was apparent that this is not consistent with our mechanistic understanding of the recharge processes.In the simple model, bypass flow was responsible for sustaining higher water levels during the summer months, when soil moisture deficits prevented any drainage from the soil.However, previous work (Price et al., 2000;Ireson et al., 2009b) has suggested that slow drainage from the unsaturated zone of the Chalk is continuous throughout the year, and sustains recharge, even during drought conditions.Thus, slow drainage, not bypass flow, sustains summer water levels.Moreover, the nature of bypass flow in the Chalk unsaturated zone has been explored (e.g. Lee et al., 2006;Ireson et al., 2009aIreson et al., , 2012;;Ireson and Butler, 2011, to name only some of the more recent studies).In a three year record at the study site used in this paper, Ireson and Butler (2011) observed 18 out of a total of 536 rainfall events led to a recognisable bypass recharge response at the East Ilsley monitoring well.Only rainfall over a threshold in volume and intensity, also dependent on antecedent soil moisture, led to bypass recharge.By contrast, in the modeling exercise in this paper, we found that a rainfall intensity threshold did not improve the model performance, and thus every single rainfall event results in bypass recharge.On the basis of these two significant inconsistencies, we conclude that the simple recharge models used widely in operational groundwater models are providing the right answers (especially those run on a monthly timestep), but for the wrong reasons.This explains why conventional recharge models for the Chalk that apply a fixed bypass fraction have Introduction

Conclusions References
Tables Figures

Back Close
Full  (1970).Our study suggests a lower value may be better (we found approx 5 % to be optimal on a monthly timestep), though this is likely to be a function of the site specific properties of the Chalk matrix and fractures.Therefore, for assessment of water resources in an average year, the existing models may be fit for purpose.However, these models may not be able to adequate predict the water table response under extreme wet or dry conditions, or under a changed climate.High intensity rainfall events have lead to significant summer recharge responses (Lee et al., 2006;Ireson and Butler, 2011) with the potential to cause rapid flooding, and to drive contaminants into the aquifers.Such responses could not be simulated with the existing models.Ireson and Butler (2011) have also shown that Richards' equation models based on a dual permeability, equivalent continuum approach are unsuitable for such conditions.Thus, modelling such responses remains a significant challenge.Above average rainfall sustained over multiple months has lead to groundwater flooding (Hughes et al., 2011).
Under such conditions we might expect the buffering capacity of the soil/unsaturated zone to be underestimated by the simple model, since with no soil moisture deficit all effective rainfall is passed directly to the water table as recharge.Whether this is a significant problem remains to be seen.Finally, during extreme drought conditions, we would expect the discrepancy between slow drainage of the unsaturated zone and bypass flow to become exaggerated.Introduction

Conclusions References
Tables Figures

Back Close
Full  Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Brouy ère Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the flux into each saturated (i.e.ψ > 0) cell (m d −1 ), z b is the elevation of the base of the aquifer (m), h is the elevation of the water table (m) which is a function of time, A is the cross-sectional area per unit width Introduction Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | tential evaporation is reduced to actual evaporation depending on soil water stress distributed over a rooting zone based on the Feddes et al. (1976) model described above.Multiple linear regression on log transformed variables identified an equation Discussion Paper | Discussion Paper | Discussion Paper | that could very convincingly (R 2 = 0.9998) reproduce the behaviour of the Feddes model, as shown in Fig. 6.The regression equation was log ) are solved in sequence, in a time loop from j = 2 Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the limitations on computing power at that time.Parameters were sampled from uniform random distributions as follows: 10 mm ≤ RC ≤ 2000 mm; 10 mm ≤ (PWP − RC) ≤ 2000 mm; 0 < BF ≤ 0.3; 0 < TH ≤ 30 mm d −1 Figure 8a also shows recharge and actual evaporation simulated by the simple and benchmark model, as well as soil moisture deficit, SMD, simulated with the simple model.These variables did not influence the calibration and are an additional test of the overall coherence of the simple model.SMD values typically reach around 0.5 Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | tended to work well.Typically 15 % bypass has been used, dating back to Smith et al.

Table 2 .
Simple model states, fluxes and parameters.