University of Birmingham Linking soil moisture balance and source-responsive models to estimate diffuse and preferential components of groundwater recharge

. Results are presented of a detailed study into the vadose zone and shallow water table hydrodynamics of a ﬁeld site in Shropshire, UK. A conceptual model is presented and tested using a range of numerical models, including a modiﬁed soil moisture balance model (SMBM) for estimating groundwater recharge in the presence of both diffuse and preferential ﬂow components. Tensiometry reveals that the loamy sand topsoil wets up via preferential ﬂow and subsequent redistribution of moisture into the soil matrix. Recharge does not occur until near-positive pressures are achieved at the top of the sandy glacioﬂuvial outwash material that underlies the topsoil, about 1 m above the water table. Once this occurs, very rapid water table rises follow. This threshold behaviour is attributed to the vertical discontinuity in preferential ﬂow pathways due to seasonal ploughing of the topsoil and to a lower permeability plough/iron pan restricting matrix ﬂow between the topsoil and the lower outwash deposits. Although the wetting process in the topsoil is complex, a SMBM is shown to be effective in predicting the initiation of preferential ﬂow from the base of the topsoil into the lower outwash horizon. The rapidity of the response at the water table and a water table rise during the summer period while ﬂow gradients in the unsaturated proﬁle were upward suggest that preferential ﬂow is also occurring within the outwash deposits below the topsoil. A variation of the is shown to reproduce the observed water table dynamics well in the lower outwash horizon when linked to a SMBM that quantiﬁes the potential recharge from the topsoil. The results reveal new insights into preferential ﬂow processes in cul-tivated soils and provide a useful and practical approach to accounting for preferential ﬂow in studies of groundwater recharge estimation.


Introduction
Many aquifer recharge estimation methods have been developed for application to water-supply, water-quality, agricultural, and ecohydrologic problems (Scanlon et al., 2002;Nimmo et al., 2005). Of these there is no single technique that can serve as a standard; each has apparent deficiencies as well as particular advantages. It is widely recommended practice to apply multiple methods to any problem requiring recharge estimation (Lerner et al., 1990;Scanlon et al., 2002;Healy, 2010). Method comparison can be useful to increase confidence in the estimates and also to highlight different features of the recharge. For example the steady state Darcian method (Nimmo et al., 1994) indicates the steady-flow component of recharge, whereas the water table fluctuation method  is sensitive to episodic and sometimes seasonally varying recharge. The augmentation of a traditional method with one that is specific to alternative flow modes could have great value. In this work we apply the multiple-method guideline for recharge estimation by combining two largely complementary methods in order to obtain an estimate that more fully represents the diverse processes contributing to aquifer recharge.
A soil moisture balance model (SMBM) is a valuable approach to recharge estimation, relying on such hydrologic inputs as infiltration and evapotranspiration (Alley, 1984). Given that water is conserved in the land-atmosphere system, recharge equals the difference between the amount of water input and the amount that goes to fates other than recharge. Significant components of the water balance typically include precipitation, evapotranspiration, soil-moisture storage, and recharge. Water transfer between SMBM components is governed by traditional principles of diffuse flow. Because recharge is computed as a residual and is frequently smaller in magnitude than the components that are differenced to compute it, all of those components must be well estimated to keep the computed recharge uncertainty from being very large.
Preferential flow in the unsaturated zone presents a difficulty for SMBMs as well as other recharge estimation methods. A common natural process, essential to plant health and ecosystem function, preferential flow bypasses much of the unsaturated medium along paths such as wormholes, fractures, and fingers of enhanced wetness, especially when water is copiously supplied. Because it typically generates highspeed, high-volume flow with minimal exposure to solid earth materials, preferential flow is a dominant influence in many problems of recharge, as well as infiltration, contaminant transport, and ecohydrology. It can produce recharge in the absence of a recognizable wetting front. It can be difficult to account for as it does not fit readily into traditional unsaturated zone flow theory, which emphasizes diffuse and equilibrium modes of flow. The key problem for a SMBM is that preferential flow can move substantial amounts of water through the subsurface in ways not amenable to a strict component-based accounting for volume of water (Ireson and Butler, 2011). For example, preferential flow may transfer newly infiltrated water directly to recharge within a time span so short that it does not allow water to first reduce soil-moisture deficits or become evapotranspiration (ET) (Cuthbert and Tindimugaya, 2010). Because standard SMBMs used for groundwater recharge estimation do not allow for processes whereby recharge can occur in the presence of a soil moisture deficit, the resulting model estimates would be in error for recharge affected by preferential flow.
Preferential flow is associated strongly with certain features of the subsurface. Its initiation can depend on threshold effects, perhaps related to input fluxes or antecedent soil moisture (Shipitalo and Edwards, 1996;Hardie et al., 2011). It can depend on small-scale features of the soil, such as macropores, whose conductance and connectivity are unknown. For example, Rosenbom et al. (2008) observed markedly different behaviour in continuous vs. dead-end macropores: during rainfall, wormholes that ran continuously down to a subsurface absorptive layer allowed freeflowing water, whereas those that terminated at shallower depths remained plugged without significant flow. The position and effectiveness of subsurface impeding layers can similarly exert strong influence. Beven and Germann (1982) noted that ploughing may cause a truncation of macropores at a particular depth. Andreini and Steenhuis (1990) found that preferential macropore flow in no-till soil flowed through the soil profile without interruption, whereas in tilled soil in which a plough pan had formed such flow occurred only below the plough pan. Direct experimental evidence (Allaire-Leung et al., 2000;Su et al., 2003) shows how it is possible for preferential flow to commence at a depth within the subsurface where a layer containing macropores underlies a layer which has none. A modified SMBM that can account for processes like these would have expanded versatility and reliability.
The preferential flow model of Nimmo (2010) is called source responsive because it allows for water at depth to respond sensitively to changing conditions at the source of water input. To avoid the need for the large numbers of parameters in most existing preferential flow models, which are often impossible to estimate a priori, the source-responsive model employs empirical relationships with a basis in laminar flow theory and properties of typical earth materials. This model employs a two-domain configuration, equivalent to a medium that has relatively finely pored matrix material (domain D) with interpenetrated macropores (domain S). In the D domain flow is by diffusive processes, quantifiable by the Darcy-Buckingham law and Richards' equation, although alternatives such as a SMBM could be used instead. In the S, or source-responsive, domain water moves by preferential flow and free-surface films, and the input of water typically dominates the flow. The S domain flow formulation is based on linear equations and thus is mathematically far simpler than Richards' equation. To predict unsaturated flow the source-responsive model requires quantitative characterization of (1) internal macropore facial area as a function of depth M(z), representing the capacity for preferential flow, and (2) an active-area fraction f (z, t), indicating how much of the preferential flow capacity is active at given depth and time. The values of M and f do not depend in general on moisture state but rather on profile-scale properties of the medium and water-input conditions (e.g. rainfall rate). Independence from local moisture state allows these functions to quantify processes that proceed without immediate inter-domain equilibration. Although it was developed around the concept of flow through macropores, the sourceresponsive model can also be useful for fingered, funnelled, or other modes of preferential flow. Much as unsaturatedflow models based on straight capillary tubes can effectively represent flow through the grossly irregular pore space of soil, this model based on macropores can represent flow through geometrically different preferential flow paths. Thus its applicability, which depends on there being some form of preferential flow, does not strictly require that this be specifically macropore flow. The source-responsive model has been applied recently in a recharge application by Mirus and Nimmo (2013).
Our specific objectives in this paper are to: 1. Evaluate soil-moisture dynamics at a field site in Shropshire, using the extensive data set available, to assess the type of diffuse and preferential flow processes that control recharge and its relation to soil moisture.
2. Test whether a traditional Darcy-Buckingham approach can adequately represent the recharge processes.
3. Develop a practical and parsimonious modelling approach based on coupling of a modified SMBM with the source-responsive preferential flow (SRPF) model to quantify the processes that are not accounted for in a standard SMBM approach.
4. Show whether the new combined model can perform adequately for the encountered field conditions.

Locality
The experimental site is located in a lowland area of gently undulating terrain in Shropshire, UK ( Fig. 1) within the catchment of the Potford Brook (catchment area 22.5 km 2 ), a tributary of the River Tern (catchment area 880 km 2 ), itself a tributary of the River Severn. It is underlain predominantly by a Permo-Triassic sandstone aquifer of regional significance with a covering of variable superficial deposits comprising glacial till, glaciofluvial deposits and valley alluvium and by soils which generally reflect the nature of the underlying superficial geology. Long-term (1970Long-term ( -1999 annual average rainfall and potential evapotranspiration for the Potford Brook catchment are around 670 mm and 597 mm, respectively . The areally averaged recharge to the aquifer in the catchment is thought to lie in the range 110-127 mm a −1 , and a previous study suggests values as high as 240 mm a −1 in outcrop/outwash sand covered areas (Streetly and Shepley, 2005). There is significant uncertainty in recharge distribution and hydraulic processes at the field scale, which forms the context for the research presented here Thatcher, 2009).
The site is located adjacent to Hollycroft Farm (NGR SJ 6408 2321) and was under winter wheat at the time of developing the site (June 2004) but was cleared by hand at harvest time (August 2004) to let grass (and other wild vegetation) become established.

Site hydrogeology
The geology of the wider catchment is described in detail elsewhere (Cuthbert, 2006;Cuthbert et al., 2009). The geology of the site is known from 4 cored boreholes, from shallow augering undertaken for the installation of the tensiometers and Time Domain Reflectometery (TDR) access tubes, and indirectly from Electrical Resistivity Tomography (ERT) surveys.
A dark reddish-brown fine to medium grained sandy topsoil (loamy sand) of 0.3 to 0.5 m thickness (ca. 0.4 m at the location of the tensiometers) covers the site. Macropores are evident in the topsoil as indicated in Fig. 2. Double ring infiltration tests indicate that the saturated hydraulic conductivity of the topsoil is in the range 0.01 to 0.1 cm min −1 (Cuthbert, 2006). No overland flow or significant ponding has been observed on the site after rainfall.
The topsoil overlies glaciofluvial outwash material persisting to between 2.45 and 2.7 m b.g.l. An ERT survey, oriented approximately east-west, suggests that the site is located on the west side of a roughly north-south oriented channellike structure with the outwash thickening to the east. The outwash predominantly comprises well-sorted medium sand with variable gravel content. Below around 1 m b.g.l. it is grey-brown in colour, but above this level and below the topsoil there is a slightly cemented orange-brown horizon of variable thickness. This is a zone of illuviation where iron, and most likely aluminium, oxides have accumulated into an "iron pan" above the level of permanent saturation. Thin (< 0.1 m) clay-rich beds are also present in some locations above 1.5 m b.g.l. Falling head tests in the piezometers at the site indicate the saturated horizontal hydraulic conductivity of the lower 0.4 m of the outwash deposits is in the range 0.1 (at piezometer BH6) to 3.8 m d −1 (Cuthbert, 2006).
Below the glacial outwash materials lies laminated brown glaciolacustrine clay persisting to 6.7 m b.g.l. on the east of the site but thinning to the west. The type of deposit underlying this clay varies across the site comprising layered till and glaciolacustrine sand in different combinations. Weathered slightly clayey, fine to medium grained Permo-Triassic sandstone underlies the superficial deposits at between 7 and 8 m b.g.l.
This paper is concerned with the groundwater hydrology of the permeable uppermost soil and outwash deposits in which a perched water table is present. A small amount of vertical leakage is estimated to occur (a few mm a −1 ) through the underlying low permeability glaciolacustrine clay deposits (Cuthbert, 2006) to the sandstone aquifer below, which is unconfined in this location.

Instrumentation and monitoring
The site was highly instrumented (Cuthbert, 2006) with the following of relevance to this paper.

Tensiometers
A nest of 5 tensiometers (T1 to T5) was installed, with specification and installation techniques as described by Greswell et al. (2009) and Cuthbert et al. (2009), respectively. Thermocouples were attached to each tensiometer just above the level of the ceramic cup. The bases of the ceramic cups were placed at 33, 73, 88, 132 and 159 cm b.g.l. for T1 to T5, respectively. The tensiometers were measured manually from August 2004 to October 2004 and logged thereafter at a 5-min frequency until July 2005. Tensiometer T2 malfunctioned within 1 month of installation due to electronic failure. The other tensiometers generally performed well, although pressure effects, presumed to be due to freezing of the water column, occurred at times correlated with sub-zero ground temperatures measured by on-site thermocouples. These data have been removed and not used for analysis. During the summer period a daily cycle of pressure variations due to heating/cooling of the tensiometer water column also caused artefacts in the data, although the general pattern of pressure changes due to moisture content variations is still clearly discernible.

TDR access tube
A 157-cm-deep access tube (M2) was installed for use with a TRIME TDR probe (IMKO, Germany). Readings were taken at an average frequency of 2 weeks at 1-to 3-week intervals for 12 months beginning in August 2004. For each monitoring period, 3 TDR moisture content values were taken at 10cm depth intervals, at orientations separated by 120 • , to enable an average reading to be calculated (see Cuthbert, 2006, for more details of the monitoring strategy).

Piezometer
A shallow piezometer (BH6) was installed within a cored borehole situated approximately 4 m from the tensiometers and TDR access tube. The monitored section (2.2 to 2.6 m b.g.l.) comprised a pre-fabricated filter pack around 25-mm ID plastic casing with 0.3-mm slots and the installation was sealed above the filter pack with bentonite, the top 30 cm being filled with a cement grout, and then fitted with a cover at ground level. It was fitted with custommade electronics including a differential pressure transducer (vented to atmosphere) to enable data to be logged automatically at regular intervals (Greswell et al., 2009). Water levels were recorded manually every few days from June to October 2004 and were then automatically logged from November 2004 to September 2005 at 5-min intervals. The piezometer transducer/data logger equipment gave results in excellent agreement with manual dip readings taken sporadically through the monitored period.

Climate data
Hourly rainfall data were available from the Bowling Green climate station situated approximately 2.5 km northwest of the site. A small number of missing data were infilled using a more distant rain gauge (from Oakley Folly, approximately 15 km to the northeast) using an adjustment factor based on a linear double mass plot. Monthly MORECS (UK Meteorological Office data; Hough and Jones, 1997) potential evapotranspiration (PE) values were also available for the locality.

General observations
Daily averages of hydraulic head across the whole monitoring period are shown in Fig. 3 and indicate that heads recorded by deeper tensiometers when saturated are indistinguishable from each other (accepting a small error in the measurements) with an almost identical pattern to the piezometer. T5 always recorded positive pressure heads. A component of lateral hydraulic gradient exists within the shallow groundwater system between the tensiometer/TDR location and the piezometer location. This suggests that there is a component of lateral drainage through the outwash deposits consistent with the presence of field drains as indicated by the landowner.
An annual seasonal cycle can be seen in the water table with superimposed recharge events separated by recessions. The levelling off observed in the water table in the summer is suggestive of minimal drainage into the underlying glaciolacustrine clay as would be expected from the measured low permeability (Cuthbert, 2006). "Winter" recessions reach a floor of around 1 m b.g.l. at the location of the tensiometers, for example during mid-March when all tensiometers imply hydrostatic equilibrium in the whole profile. Steeper and deeper recessions occur in the summer period, levelling off at around 1.4 m b.g.l. T1 is much more responsive than the deeper tensiometers, and a downward gradient persists during much of the monitored period, except during the height of the summer. Notably, the piezometer and tensiometers, excluding T1, are in hydrostatic equilibrium during the whole monitored period except during the height of the summer periods.

Summer responses
PE reaches a maximum in June, averaging approximately 3 mm d −1 (Fig. 3), greatly increasing evapotranspirative demand for soil moisture during these summer periods. As a result, pressures in T1 and T3 are frequently negative, reaching tensions of a few hundred cm of water, although during prolonged/intense summer rainfall events tensions in the topsoil are severely reduced. Water table rises are seen under low  water table conditions, for example during late August 2004 (5-cm rise) and late June 2005 (8-cm rise). During the former event it is notable that the water table responds while an upward flow gradient is present in the profile as shown by a lower head recorded in T3 than in the lower tensiometers (T4 and T5).

Autumn to spring responses
A subset of high temporal resolution (5-min) data for March 2005 is shown in Fig. 4 and is illustrative of the behaviour of the system outside of the summer periods. T1 (33 cm b.g.l.) shows complex threshold behaviour responding almost immediately to certain rainfall events but showing no response to others; e.g. the first two rainfall events on 22 March produce no response, but then the third, of smaller magnitude, produces an effect within minutes in T1 with the tension almost reducing to zero. After rainfall ceases a "background" tension is quickly re-established but of lower magnitude than that preceding the event. Between such events, T1 shows some recession (e.g. between 25 and 30 March). The deeper tensiometers do not show responses to all events in which T1 responds unless tensions in T1 are reduced for long enough so that the "background" tension is reduced to close to zero (< 10 cm or so). Once this threshold is reached, a water table rise occurs within minutes (e.g. water table rise starting from around 1 m b.g.l. on 31 March) followed by a recession seen in all tensiometers. These types of responses are shown consistently by the data for all the main recharge events throughout the monitored period.

Soil moisture data
The general pattern of moisture content changes evident from the TDR data is consistent with the hydrodynamics inferred from the tensiometer data (Fig. 5). In particular, a zone of low moisture contents between 40-60 cm b.g.l. is consistent with the compacted base of plough/topsoil seen in the geological logs. However, there is large uncertainty in the absolute values in moisture content due to the lack of site-specific calibration for the instrument (Cuthbert, 2006). Given the low temporal resolution of the data, the TDR data cannot be used to infer much more information about the hydrodynamics than the pressure data in this instance.

Conceptual model
The key observations so far described lead us to form the following conceptual model of the hydraulic processes operating at the field site. Fig. 6 summarizes this visually, as well as incorporating the main features included in the numerical model described later. The profile above the low-permeability glaciolacustrine deposits (at 2.7 m b.g.l.) can be subdivided into three main horizons comprising (1) "topsoil" -loose organic rich sandy material, which is disturbed by ploughing, (2) a partially cemented and compacted "plough pan", and iron pan material of lower matrix permeability derived from (3) "outwash"unconsolidated glaciofluvial sand and gravel materials. The shallow groundwater system drains laterally to field drains, with minimal drainage vertically into the underlying clay.
It is apparent that preferential flow processes are in operation in the topsoil from the extremely quick responses observed in the tensiometer T1. T1 is most likely associated with macropore(s), which become active only if rainfall intensity and antecedent conditions are favourable. Once the macropores are active, flow occurs rapidly to the level of T1 while "abstraction" (in the sense of Hincapié and Germann, 2009) of water is also occurring from the macropores to the soil matrix. The pan horizon is likely to have a much lower vertical matrix hydraulic conductivity than the overlying topsoil. Thus persistent hydraulic gradients between the topsoil-pan layer and pan-outwash layer may not necessarily be evidence of significant matrix flow. The lack of water table rise during the last 10 days of March in the presence of strong vertical hydraulic gradients in the upper profile is evidence of this. Despite the topsoil responding to rainfall events during the rest of the month with positive pressures building up, no response in any of the deeper tensiometers is seen until the large rainfall event on 31 March which initiates a recharge event as described below.
The rooting depth of most crops grown at the field site is significantly greater than the depth of the topsoil. Hence it is reasonable to expect that old root networks may provide preferential pathways through the plough pan and upper outwash zone but also that seasonal ploughing activity is likely to truncate macropores developed in the topsoil from any established at deeper depths in the profile (Beven and Germann, 1982;Frey et al., 2012). This vertical discontinuity enables the macropores eventually to fill at the level of T1, but once the supply of water to the macropore reduces sufficiently, or ceases entirely, the redistribution of moisture from the macropores to the matrix causes a quick recession of tension back to the background level of the matrix. This background level of tension may be lower than that preceding the event due to the transfer of water from the macropore to the matrix during the event.
Furthermore, if the matrix tension of the upper soil profile above the level of the macropore truncation is reduced sufficiently (i.e. saturated conditions at the base of the topsoil) then rapid flow occurs from the topsoil to the water table. Thus the upper soil horizon is a "source" to which the water table is "responsive". Two observations suggest that flow through the outwash material may be occurring, at least partially, via preferential flow processes. First, the water table response during late August 2004 (Fig. 3) occurred while an upward flow gradient was present between T3 (85 cm b.g.l.) and T4 (129 cm b.g.l.). Second, water table rises are very rapid and occur without obvious departure from hydrostatic conditions in the unsaturated part of the outwash horizon. Once the upper layer is drained sufficiently, preferential flow in the lower layer ceases and the whole profile exhibits recession due principally to lateral drainage (the underlying glaciolacustrine clay is of very low permeability). During much of the year, apart from the warmest summer period, the water table is very shallow and the whole profile close to saturation and water transfer to the matrix during preferential flows within the outwash deposits is minimal. During very dry periods, when the evapotranspirative demand is higher and rooting depths greater, moisture may be removed directly from both the topsoil and the upper part of the outwash deposits, and an upward flow gradient may develop within the outwash horizon, contributing to the water table recession during this period. Summer storms may quickly activate macropores in the topsoil, and preferential flow pathways in the lower profile may also become active under prolonged rainfall enabling groundwater recharge to occur. In this case water transfer to the matrix of the outwash horizon may be higher depending on the deficit that has accumulated in this zone, and actual recharge may be less than the potential recharge leaving the topsoil.

Hypothesis testing: what type of flow process
dominates the recharge response?
The conceptual model described above uses several lines of evidence for preferential flow through both the topsoil and the outwash deposits beneath. We consider that the tensiometer responses in combination with the direct observations of macropores within the topsoil are conclusive in this respect for this layer. The significance of preferential flow within the underlying outwash deposits compared with matrix flow and the mechanisms that cause it are less immediately apparent and require further investigation. The summer water table response in the presence of upward gradients as evidenced by the tensiometers seems to be good evidence that preferential flow does occur at certain times unless lateral groundwater flow can occur as a mechanism for inducing water table rises under these conditions. Under higher moisture content conditions, it is not obvious whether the observed quick water table responses can be accounted for using conventional unsaturated zone theory. In this section we explore the hypothesis that water table responses can be explained using 1-D Darcy-Buckingham type flow through the outwash deposits. We have achieved this through the use of the variably saturated flow modelling code FAT3D-UNSAT to model 1-D water flow through the profile with Richards' equation using hourly forcing climate data described in Sect. 2.3. FAT3D-UNSAT solves Richards' equation for heterogeneous single porosity media using a block-centred finite difference formulation with backward differences in time and central differences in space. Soil properties are described by the van Genuchten (1980)- Mualem (1976) soil moisture characteristic equations. Boundary conditions are assumed to be piecewise constant over timesteps. A Newton-Raphson scheme with damping is employed to solve the resulting non-linear equations. Time stepping is automatic to ensure stability and convergence of the solution. The code has been extensively verified against HYDRUS (Šimůnek et al., 2011) and analytical models for a range of problems and is found to provide excellent agreement with the comparison solutions. In addition to piecewise general head and flux boundary conditions, FAT3D-UNSAT can model seepage boundary conditions, and variable evapotranspiration controlled by soil moisture content. Automatic calibration is not possible with the tool, and a process of manual calibration for the modelling was undertaken. A model grid was employed with uniform one-centimetre vertical cell resolution to model the top 2.6 m of the profile to the depth of the base of the outwash deposits. The initial and boundary conditions used for the analysis of the field site are as follows.
The initial modelling period for this test was chosen to be from 1 March to 5 April, while the calibration period was set to be from the 19 March to 5 April. The lead-in period of 19 days was selected to allow a simple initial condition to be employed and any errors in the initial conditions to be dissipated before the calibration period. A uniform variation of pressure head with depth was used with the water table at −0.96 m below ground level for the initial condition, with a gradient of −0.9 representing a gentle downward flux condition. The calibration period was insensitive to alteration of the initial condition.
Boundary conditions were applied as follows: the hourly precipitation from the nearest rain gauge was added to the top boundary node uniformly over each hour. Potential evapotranspiration was applied uniformly over each calendar month, based on the available MORECS data. Actual evapotranspiration was calculated for each node of the top 10 cm of the soil profile, based on the moisture content. For the calibration period AE = PE.
The top node was set to remove excess pressure head by lateral drainage when the pore pressure exceeds zero. In the saturated zone, lateral flows were implemented using a general head boundary condition for all nodes from 1.2 m bgl to the base of the outwash. Physical properties are defined for a range of soil types, and these can be assigned to the cells of the model to create a layered model.
A variety of layer thicknesses and parameter combinations were used to attempt to simulate the observed response in the tensiometers and piezometer. The initial choice of the van Genuchten characteristics were made using simple comparisons between the particle size distributions of the soils at the field site and the USDA standard soil classes. These were then varied to yield the final fits. The use of four layers for this exercise reflected the different visual descriptions for the upper and lower glacial outwash deposits used to guide the initial choice of van Genuchten parameters. Based on a careful sequence of calibrations it could be shown that there are many layering and parameter combinations that give similar responses, but none were found that simulated well the  Fig. 6. Summary of conceptual and numerical processes for the SMB-SRPF model. In the flow diagram, the left side with blue arrows indicates the source-responsive processes that supplement the traditionally recognized moisture balance processes on the right side.
responses seen in T1. A typical "best fit" is shown in Fig. 7a. The inability for a relatively complex single domain 1-D Richards' equation model to simulate the observed hydraulic response in T1 supports our qualitative interpretation that preferential flow processes are operating in the topsoil. Better fits were obtained for the lower tensiometers (T3-T5) and piezometer (BH6) in the outwash deposits. In particular, the model is able to produce sharp water table rises such as that shown in Fig. 7a, as long as the retention capacity of the outwash is sufficiently high (a high value for the van Genuchten n parameter). However, for models that achieve this feature of the observed data, "unwanted" water table responses are simulated during other periods, for example the significant early rise of the water table beginning on 24 March 2005 in Fig. 7a. When the model developed for the one month was tested for the full year (July 2004-June 2005), the results were very poor and a recalibration exercise was undertaken for the full year. For this recalibration the effort was concentrated on returning the observed water table response. Figure 7b presents the results of the recalibration. It proved possible to reproduce the water table variations with reasonable accuracy, but at the expense of reproducing the observed conditions in the topsoil layer.
Two observations can be made based on the results of this modelling. First, the results for the summer periods produce the main water table rise events and in these periods water tables do rise in the presence of upward gradients in the overlying unsaturated zone. This is because the RE model invokes a general head boundary, which allows water to move both towards the monitored site as well as away from it. Essentially the model is allowing for the possibility of persistent non-horizontal groundwater gradients in the region of the test site. Second, it has been possible to match the winter period well for either the interval (October to January) or the interval (February to April) but not both. Moreover, it has not been possible to retain the sharpness of the transition from rising to falling water table levels observed in the data in the wet period. The model does recreate elements of the rising water levels at the end of January 2005 even in the absence of heavy rainfall. Thus we cannot completely reject the hypothesis that water table responses can be explained by solely invoking 1-D Darcy-Buckingham type flow through the outwash deposits, but the results do suggest that preferential flow below the topsoil is contributing to the recharge responses observed over the seasons. This modelling exercise also proved useful in demonstrating the significance of the pan layer discontinuity in restricting/moderating flow from the overlying morepermeable topsoil, consistent with the conceptual model described above.

The need for simple preferential flow models
Based on the discussion so far, it is clear that a complex set of hydraulic processes control recharge to the shallow water table at the field site. It is likely that both capillarydominated matrix flow and preferential flow processes interact to give the observed response. Modelling these processes using state-of-the-art dual permeability models requires many input flow parameters for each medium type.  and van Genuchten (1993) approach requires 16 parameters (or possibly 11 if various other assumptions are made), or if a kinematic wave approximation is made for the fracture domain then 10 parameters is sufficient. If no flow is assumed in the matrix domain (i.e. a "dual porosity" formulation), between 9 and 11 parameters are needed for each soil type. Thus, to model the profile described in this paper would require between 18 and 32 parameters for a two-layer model (e.g. topsoil and outwash) or 27 to 48 parameters for a 3-layer case (e.g. topsoil, plough pan, outwash), just to characterize the water flow (additional parameters are needed for other aspects such as root water uptake). Such models and their inherent computational demands are of little use for estimating groundwater recharge at the field and catchment scale, and more parsimonious approaches are very desirable. We propose the following approach.

Model structure
The model couples a modification of an established soil moisture balance model of the topsoil with a simplified preferential flow model of the outwash deposits as a two-layer system. A schematic outline is shown in Fig. 6.
The SMBM for the topsoil is modified from that developed for use in regional groundwater models in the UK, as described by Hulme et al. (2001) and Rushton et al. (2006), which makes use of a combined crop co-efficient approach (K c ) and other concepts taken from Allen et al. (1998). The total available water (TAW) in the soil is defined as (1) where θ FC and θ WP are fractional soil moisture contents at field capacity (FC) and wilting point (WP), Z r is the rooting depth of crop, Z e is the thickness of the soil layer subject to drying by evaporation, B = fractional area of bare soil (i.e. crop absent). The readily available water (RAW) is defined as where p is a factor normally between 0.2 and 0.7 (Allen et al., 1998;Table 22). If the soil moisture deficit in the topsoil (SMD s ) is greater than the RAW, then the actual evaporation rate (AE) is reduced using a stress co-efficient (K s ) as follows: where potential evapotranspiration PE = K c PE 0 and PE 0 is the reference crop (grass) potential evapotranspiration rate. Using an input time series for rainfall (RF) and PE 0 , the model algorithms calculate time series of AE and the rate of potential groundwater recharge (R pot ). Overland flow has been assumed to be zero with all rainfall becoming infiltration. On days where the soil is under stress and rainfall occurs that is less than PE, the rainfall is transpired plus a further amount from the soil equal to the remaining evaporative demand modified by the stress co-efficient. If rainfall exceeds the PE, then the excess reduces the SMD. If the SMD becomes negative, this water is accumulated in a store (of amount D) and released as potential recharge (R pot ) to the outwash deposits according to a limiting flux based on Nimmo's (2010) source-responsive model: where V u L u is the product of the film flow velocity and film thickness, M lim is the macropore facial area density and t is the model timestep. A constant default value of V u L u = 5.5 × 10 −10 m 2 s −1 was assumed, with M lim treated as the controlling variable.
If, after satisfying the demands of the topsoil, there is a remaining evapotranspirative demand (PE r ), this is taken as capillary rise (CR) from the water table (at depth b.g.l. = GWL) within the outwash deposits using the following algorithm based on an extinction depth (L) concept: Any evapotranspirative demand still remaining is accumulated as a soil moisture deficit (SMD o ) in the plough pan layer up to maximum amount SMD oLim . While SMD o is greater than zero, the rate of actual recharge (R act ), is calculated as follows: where A is a constant factor between 0 and 1, or else R act is simply equal to R pot . This can be thought of as a "bypass" mechanism whereby a constant proportion of the potential recharge available from the topsoil bypasses the plough pan layer, becoming actual recharge. Alternatively this could be considered as a very simple dual porosity formulation for governing preferential flow whereby an amount equal to (1 − A) · R pot is transferred from preferential flow pathways to matrix contingent on the presence of a moisture deficit in the matrix. It is assumed that the vertical loss to the underlying clay is insignificant and thus the elevation of the water table (WT) above the elevation of the field drains (at depth = FD) is controlled by the balance of lateral drainage, capillary rise and groundwater recharge as follows: where t is time, τ is a recession constant and S y is the specific yield, and WT = FD − GWL. The resulting number of hydraulic parameters controlling the recharge behaviour (i.e. ignoring the crop aspects) is just six (θ FC , θ WP , M lim , L, A, SMD oLim ) with an additional 3 parameters for controlling the resulting water table fluctuation (τ , S y , FD).

Model results and discussion
The model was run on an hourly timestep using forcing climate data for exactly one year (1 July 2004 to 30 June 2005) as described in Sect. 2.3 with initial parameter estimates and then refined against the observed groundwater levels. The model refinement procedure was as follows: 1. During the winter period (November 2004 to April 2005) AE was always equal to PE and the water table responses were effectively governed by the forcing data and just three parameters -S y , τ and M lim . The model was run for this period with an initial proportion of bare soil 0 to 1 depending on cover -SMD oLim maximum soil moisture deficit in outwash 28 mm L extinction depth for capillary rise 1.5 m A preferential flow abstraction factor 0.5 -FD depth of field drains 1.1 m b.g.l. GWL i initial groundwater level 1.2 m b.g.l. SMD i initial topsoil moisture deficit TAW mm SMD oi initial outwash moisture deficit SMD oLim mm water level equal to the observed water table elevation and, as a starting point, optimized to maximize the Nash-Sutcliffe efficiency. As is often the case with this type of optimization, multiple parameter sets can yield similarly good fits by this criteria, but it was found that the model fit was always poor during February 2005. Assuming that the specific yield of the outwash materials does not change with time, it is physically impossible for rainfall as measured at Bowling Green to generate the observed water table rises at the site even in the absence of any evapotranspiration during this period. This is in contrast to the Richards equation model, which does produce a better fit for this period. This illustrates how the different model forms may characterize some but not all of the features in the actual data. Removing this period from the optimization enabled a fit to be made with a higher Nash-Sutcliffe efficiency of 0.9 for a physically realistic value of S y as given in Table 1. 2. Having fixed these parameters, the whole period was then modelled. Soil/crop parameters (θ FC , θ WP , Z e , K c , p, SMD oLim , B) were estimated using mid-range values directly from Allen et al. (1998), and the remaining unconstrained parameters (L, A) varied manually to achieve a best fit.
The best fit is shown in Fig. 8 with corresponding model parameters given in Table 1 having a good Nash-Sutcliffe efficiency of 0.89. Despite the complex processes occurring in the topsoil, the SMBM approach is very effective at predicting the timing of the main recharge events. The magnitude of the recharge events is as uncertain as the forcing meteorological data, but the fitted value of around 6 % for S y is certainly physically plausible given the nature of the outwash materials and is consistent with the Richards equation model. The 676 mm of rainfall modelled over the year is distributed by the model into 559 mm of AE and 123 mm of actual recharge. The PE and potential recharge for the period were 624 mm and 172 mm, respectively. Lateral flow from the outwash deposits was 91 mm, capillary rise was 34 mm, and the differences in the various inflows and outflows corresponded to storage changes in the soil and outwash deposits. The general shape of the groundwater recessions are also simulated well using the combination of a recession constant for lateral flows and extinction-depth-controlled capillary rise for vertical losses. Furthermore, given the extremely simple representation of exchange between the preferential and matrix domains used, the model is able to simulate reasonably well the summer recharge events which occur in both years using a value of 0.5 for the bypass factor A. During the winter, the shape of the water table rise is most sensitive to the parameter M lim , which controls the rate of drainage of potential recharge from the topsoil into the outwash deposits. Figure 9 shows the detailed response for two large recharge events simulated by the model, and sensitivity to this parameter. The model works well for M lim between 250 and 750 m −1 . There is a suggestion from the data that larger recharge events lead to steeper water table responses. With a longer data series to work from, it may be possible to re-structure the model to accommodate a variable M value and also introduce a time-varying active area fraction as suggested by Nimmo (2010). As in the water table fluctuation case study presented by Nimmo (2010), our model's assumption of immediate water table response to water arriving at land surface neglects the time lag inherent from the finite speed of transit through the unsaturated zone. For a response dominated by preferential flow, however, the assumed zero transit time may be an adequate approximation for a response time that is orders of magnitude faster than the other, nonpreferential, flow processes occurring through the medium. For greater realism, a finite time lag could be incorporated to correct for this, though at the cost of increased model complexity. We consider that the model as presented here performs acceptably given the simple parameterization implemented.
In order to illustrate the merits of the new model and the differences among the three models, we have compared the recharge results. While the overall estimation of recharge among the models only varies by around 8 % (SMB-SRPF -123 mm; RE -126 mm; SMBM -117 mm), substantial differences in behaviour are seen, particularly during wet periods. Comparative results shown in Fig. 10 indicate the distinctively different behaviour of each model over the year. Both the SMB-SRPF and the RE models reproduce the dry-period recharge events during later summer 2004. The SMB-SRPF model shows point recharge responses during the winter, while the RE model presents a more continuous recharge response. Recharge rates are consequently much lower than those for the SMB-SRPF model. The peakiness of the SMBM-SRPF GWL response compares to the smoother RE GWL response and suggests that the results lie between these two extreme models. The RE model relies in part on lateral flow contributions to the groundwater dynamics to obtain the rising groundwater level responses in the dry period. This mechanism cannot be excluded given the available data and highlights the need for lateral monitoring at any site where horizontal flow components cannot be discounted from the water balance model. Based on the results of all models, the adoption of a basic SMBM to characterize the top soil appears to be well justified.
Other approaches to modify SMBMs to account for preferential flow have been suggested in the literature. In particular, the use of a one-layer SMBM with an additional bypass recharge mechanism is standard practice for implementing within-UK regional groundwater models, for example using the code 4R (Heathcote et al., 2003). The concept of bypass recharge has a long history beginning at least with Mander and Greenfield (1978) (cited in Butler et al., 2012) who used a constant proportion of 15 % of rainfall to bypass the soil zone to account for groundwater level fluctuations in the summer months. Typically, a rainfall threshold is used above which a constant proportion of rainfall may become recharge. Such an approach was tested for this site and found to be rather unsuccessful. Although summer water table responses could be simulated to some extent using a rainfall threshold for bypass flow, the end of summer rise and the beginning of summer recession could not simultaneously be modelled well. A greater sensitivity to antecedent moisture conditions is required than this type of simple bypass model accounts for. The new model developed in this paper represents a significant step forward by putting the modelled bypass concept on a more physical basis, informed by the observed field data.

Summary and conclusions
Examination of recharge-related processes using previously unpublished data at a field site in Shropshire, UK, reveals significant influence of both diffuse and preferential flow. Much of the time, hydraulic gradients and water fluxes at this site follow expected patterns based on traditional unsaturated flow theory. Sometimes, however, the data show behaviour not explainable by this theory alone. For example, a rising water table that indicates positive recharge sometimes occurs, while measured hydraulic gradients in the unsaturated zone suggest flow should be upward without the adoption of a lateral flow source. A relatively complex singledomain 1-D Richards' equation model is incapable of consistently simulating the observed hydraulic response in the upper horizon over a period that includes evidence of preferential and diffuse flow, thus supporting our qualitative interpretation that preferential flow processes are operating in the topsoil. Alteration or supplementation of traditional unsaturated flow models may be needed even where there is no direct observation of preferential flow. The time and space resolution of the Shropshire field measurements, even though they are much better than what is usually available, still fall short of being able to capture the effects of preferential flow, which can cause recharge at a metre or more in depth in less than one hour. The paths of such flow may be narrow and relatively fast flowing so as not to cause a sensible increase in water content in some of the layers the flow passes through.
Reasonable fits with the Richards equation-based model were obtained for certain tensiometer and piezometer data, and for sharp water table rises, but the Richards model generated unwanted water table responses during other periods. It was essential with this model to use unrealistic parameters for the top soil to gain the correct water table responses. Much more realistic values for the soil moisture characteristics could be introduced for the underlying glacial outwash when the model was calibrated to the full year of groundwater level data. The rapid transitions between pressure rise and fall in the deeper profile could only be partially reproduced with this model. Even with multiple parameters optimized for best fit, the Richards equation-based model does not give a consistent representation of the behaviour of the unsaturated flow system.
A more complex, multi-porosity model based on Richards' equation may improve the simulated response, but the benefits of this for recharge estimation for the site appear to be small. Moreover, the additional complexity would require significantly greater data collection to obtain the required model parameter values. The comparison with the simplified models suggests that there may be features of the recharge that are not entirely captured by the simplified models, and further work may be needed to establish the significance of these, but in terms of recharge magnitude and general timing the simplified models are applicable to the site conditions. Soil moisture balance models are attractive for estimation of recharge, being straightforward to apply and having a sound basis in the conservation of water within the land-atmosphere system. But neglecting certain important processes, such as preferential flow, a SMBM may often require some augmentation or refinement to produce reliable recharge estimates. We developed a practical and parsimonious modelling approach for aquifer-recharge estimation based on well-known SMBMs coupled with the recently developed source-responsive model for preferential flow. The source-responsive preferential flow model is in several ways an advantageous choice for augmenting the SMBM approach. Like the SMBM, the SRPF model is easy to implement, parameterize, and compute. It attends directly to factors that play a strong role in preferential flow, like the temporal character of water input. It is conceptually compatible with SMBMs in recognizing that the presence or absence of a source of water in various positions within the soil profile dictates the hydrologically important flow phenomena, often more than do the traditional Darcian driving forces and hydraulic characterizations of the medium. The combined SMB-SRPF model accounts for aspects of both soil architecture and preferential flow. Furthermore it is parsimonious, having just six parameters controlling water flow (θ FC , θ WP , M lim , L, A, and SMD oLim , excluding the crop parameters).
Despite the complexity of active hydraulic processes at locations such as the Shropshire site, the SMB-SRPF model is effective at predicting the timing of the main recharge events. The general shape of the groundwater recessions are also simulated well using the combination of a recession constant for lateral flows and extinction-depth-controlled capillary rise for vertical losses. With an extremely simple representation of exchange between the preferential and matrix domains, the model simulates the observed water table fluctuations well, even in the summer. The improvement in model results with explicit incorporation of preferentialflow processes further demonstrates the importance of the hydraulic discontinuity between the topsoil and underlying deposits created due to ploughing. The existence of many similar physical systems suggests that the model may have wide applicability.
This combination of models also holds promise for improvements in the SRPFM of Nimmo (2010). In particular it should help to explain the effects of thresholds and antecedent water on preferential flow, which are handled directly in SMBMs but not in the SRPFM. The successful combination demonstrated here also shows how layers within the subsurface can serve as sources for water supplying preferential flow, much like the surface-applied water sources considered by Nimmo (2007).
Further work will refine and extend the approach to other field applications and allow for upscaling to improve estimates of recharge timing and magnitude at the catchment scale. Thus insights and concepts from this combined modelling can also lead to progress in the treatment of general preferential flow and water resource issues, in addition to more realistic accounting for the diverse modes of unsaturated flow affecting recharge