Introduction

Dune system hydrology

Dunes are an important feature of the British coastline. In England and Wales, they often form shallow rain-fed aquifers draining to the periphery of the dune system; to the fore dunes and the beach and inland to low-lying land (Fig. 1).

Fig. 1
figure 1

Generic cross section through a dune system showing main groundwater and surface water features (Stratford et al. 2013)

Full size image

They typically consist of wind-blown sands (from exposed sea bed in the last glacial period) that accumulated on top of impermeable estuarine or glacial clays, which in turn form the base of the dune system aquifer and provide a distinct hydraulic barrier to the underlying bedrock strata (Robins et al. 2013; Stratford et al. 2013). The British climate provides an average annual excess of rainfall over evaporation during winter months which causes an accumulation of water within the dune systems, creating a so called “groundwater dome”, although the outflow conditions at the edges vary. The extent and height of the dome, and hence the groundwater head distribution within the dune system, is controlled by various factors including recharge, aquifer geometry and the physical properties of the sand aquifer, specifically by dune sand storage coefficient (S) and the distribution of permeability or hydraulic conductivity (K). S and K are often approximated from grain size analysis (Robins 2013) as conventional pump tests methods cannot easily be applied in thin, high conductivity aquifers where the transmissivity (defined as the product of aquifer thickness and K) declines rapidly as the drawdown from pumping increases. Moisture content profiling can be applied to help improve understanding, by providing a recharge value against which formation parameters such as K and S can then be validated (Stratford and Palisse this issue).

Geometry and spatial dimensions of the dune system are also important in controlling groundwater head distributions. The elevation of the aquifer base in relation to main drainage features, for example, controls the saturated thickness of the aquifer (b) and hence its transmissivity (T). T is defined as the rate of flow through a unit width of aquifer under a unit hydraulic gradient and is also related to the aquifer’s hydraulic conductivity (T = b*K). Other important factors controlling groundwater head distributions within dunes are the functionality of man-made drainage systems and the effective precipitation (i.e. the amount of water that infiltrates into the dunes following interception of rainfall on the plant canopy, recharge of any soil moisture deficit, evaporation from wet plant surfaces and transpiration by vegetation).

Each dune system has its own unique drainage characteristics (e.g. direct flow to the sea, free drainage inland, agricultural drains, controlled water levels or adjacent river with varying flows). Discharge to the beach beneath the fore dune creates a diffuse brackish zone rather than a saline wedge. This zone is flushed twice daily by the sea tides and, at low tide, a mixture of sea water and dune groundwater can often be seen discharging to the foreshore.

Land use varies between the systems, some are actively grazed whereas others are left un-managed or have partial tree cover. The sensitivities of the water tables may differ; one site may be susceptible to changes in vegetation, another to coastal erosion and accretion or to the rainfall regime. It is important to understand these sensitivities in order that each coastal dune system can be managed to best advantage.

Dune system ecology

Coastal dunes are valuable ecological reserves within which a series of hydrological and hydroecological environments exist (Davy et al. 2006; Davy et al. 2010; Lammerts et al. 2001). The dune habitats are strongly influenced by hydrology and groundwater levels within the dune systems especially in the ‘slacks’, the low-lying areas between dune ridges. The water table in the slacks is near the surface, or intercepts it during periods of high water levels, but has a high seasonal variation’ (Grootjans et al. 1998). Groundwater response to recharge is different in different slack types (Lammerts et al. 2001) and lags between rainfall and groundwater level change in some systems can be 4–7 months (Jones et al. 2006). Groundwater table variations determine the abundance of many slack species, in particular in wet (humid) slacks (Willis et al. 1959a; Willis et al. 1959b). Seasonality is also important, with spring water levels controlling the breeding success of European Priority species. For example, Natterjack toads Epidalea calamita, which are on the ICUN red list of threatened species (http://www.iucnredlist.org/details/54598/0), need shallow freshwater ponds in sandy environments to breeds successfully. The ponds need to 30-100 cm deep between the months of April to June, and ideally should dry out as the summer progresses (Denton et al. 1997). Should the water table in these dune systems recede, these slacks will evolve from humid to dry accompanied by a corresponding deterioration in habitat variety and quality (Ranwell 1959). In order to preserve the preferred humid slack environments, dune managers need to have some understanding of likely changes in the elevation of the water table and the factors that control them.

Groundwater monitoring in coastal dune systems

Dunes have been regularly monitored (in some cases for over 40 years) using dip wells which are usually recorded on a monthly basis (Clarke and Sanitwong Na Ayuttaya 2010) although more frequent water level monitoring (e.g., at 30 min intervals) is being undertaken at some UK dune systems, e.g., at Formby Point and at Ainsdale. These data provide information on seasonal fluctuations and inter-annual changes in groundwater levels and also the short-term response of groundwater to differing patterns of rainfall and evaporation. No clear climate change signal has yet been detected in these observations, yet climatic change and the impact of land management is a major concern for nature reserve managers and other environmental organisations (Provoost et al. 2011).

Groundwater monitoring within coastal dune systems is driven by the need to provide appropriate data and understanding for effective dune system management. Spatial and temporal data are needed to develop a conceptual model of the hydrological/ hydrogeological functioning of the dune system. This conceptualisation can in turn be tested by numerical modelling once sufficient data are available for this purpose. The aim of any groundwater monitoring programme is, therefore, to arrive at an understanding of how a dune system operates, the changes that have occurred over a range of timescales and how those changes have impacted the hydrology and ecology of the dunes system. This knowledge informs how changes in environmental drivers or dune management may impact the system in the future.

An important aspect of groundwater data collection and monitoring within coastal dune systems is to understand and quantify the different components of the hydrological water balance (Fig. 1). However, there are differences of viewpoint on the value of groundwater monitoring and modelling between the scientific community and the organisations who own and manage the dune systems. Hydrologists are interested in developing process models to assess topics such as long term climatic change effects and the effects of sea level rise whereas a managers usually want to understand the causes of shorter term inter annual variability of aquifer recharge and the impacts of changing land use policy such as grazing or removing tree cover. Nevertheless, an important aspect of groundwater data collection and monitoring is to understand and quantify the different components of the hydrological water balance (Fig. 1). This includes understanding the temporal and spatial distribution of groundwater heads (such as water table and pond levels), flow directions, rainfall recharge inputs as well as aquifer discharge, i.e. inland seepage or groundwater flow towards the sea.

Groundwater modelling of coastal dune systems

While water level observations provide valuable information, they alone do not explain the processes that drive changes in groundwater levels. A key way of dealing with this problem is the development of dynamic groundwater models, either conceptual or distributed numerical groundwater flow models. Both are data intensive techniques, but a validated groundwater model would enable ‘what if’ scenarios to be tested and the model to be tested on likely impacts of human or natural changes.

A number of different modelling approaches of differing complexities have been applied to understanding groundwater flow and solute transport in dune systems. The simplest is a water balance approach in which inflows and outflows are calculated and compared to storage changes in the aquifer (e.g. Jacobson and Schuett 1984). Analytical solutions, exact solutions of governing equations, have been applied at some sites and examples include a comparison of deterministic models with analytical element approaches which have been undertaken to simulate abstraction for the Amsterdam Water Supply (Olsthoorn 1999).

Simple regression or other “black box” approaches, in which processes are not explicitly represented, can be applied but must be used within observed limits (e.g. Wheeler et al., 2015). Otherwise, they may predict water levels that are physically impossible by extrapolating a fitted line between rainfall and water table levels. Another approach might be to develop a point or lumped model (e.g. Mackay et al., 2014) where the whole system is considered to be uniform, i.e. with constant vegetation cover and uniform distributions of grain size and permeability. This has the advantage of simplicity but a disadvantage is that dune systems are heterogeneous. A more usual approach is to develop a model that is spatially distributed (e.g. MODFLOW model of St. Fergus dunes; Malcolm and Soulsby 2000). The disadvantage is that gaps in data need to be identified and addressed, although the modelling process itself can identify these gaps (e.g. Jackson et al., 2005). Understanding spatially distributed parameters will lend itself to a regional distributed groundwater model.

A specific challenge in dunes system modelling is the representation of slacks (and groundwater-fed ponds) which can form a direct outflow from the groundwater system through evaporation (Winter and Rosenberry 1995) and, hence, can significantly impact the overall water balance. They are by their nature ephemeral, appearing in winter after periods of recharge and disappearing during summer when groundwater heads are low. As such they are similar to groundwater fed, ephemeral lakes and can be modelled using approaches for modelling Eskers (Ala-aho et al. 2015), groundwater-fed lakes (Lubczynski and Gurwin 2005) or Turloughs (Gill et al. 2013). However, the challenge for modellers is dealing with groundwater when it appears above ground level, particularly translating it into a surface water system (pond or lake), groundwater inflow and accounting for evaporation. Potential approaches include using fully-featured models such as Mike-SHE (e.g.(House et al. 2015), MODFLOW lake package (Merritt and Konikow 2000) or combining bespoke pond models with groundwater models using dynamic model linking techniques such as OpenMI (Knapen et al. 2013).

Given the proximity to the sea salinity can be an issue in terms of groundwater flow and abstraction in relation to saline intrusion (Essink 2001). A number of studies have been undertaken using variable density modelling to address this, particularly where groundwater is being abstracted from the dune system. The coastal Belgium dune systems have been extensively researched and a number of models created to help understand the challenges of sustainably extracting groundwater from dune systems. Vandenbohede et al. (2009) report the use of MOCDENS3D to simulate how groundwater abstraction in Westhoek, Belgium affects saline intrusion from the sea into a dune groundwater system.

The aim of this paper is to demonstrate how differing approaches to developing groundwater models can lead to a) better understanding and improved conceptual models of the groundwater systems, b) identifying the key drivers that affect water table levels and changes such as recharge and boundary conditions, c) the suitability of existing data to validate and test the conceptual models, including identification of additional data requirements and d) being able to simulate adequately seasonal, inter annual and longer term trends in groundwater levels.

Groundwater level observations have been recorded at over 10 dune systems in the UK and Northern Ireland, including Sefton, Braunton, Newbrough, Sandwich, Seascale, Whiteford, Merthyr Mawr, Point of Ayr, Magilligan and Portstewart. We selected two sites for our case studies which have the longest continuously monitored records of groundwater levels in the UK – Braunton Burrows in Devon, South East England (40 years) and at Ainsdale in Merseyside, North West England (44 years). Models for the two groundwater systems were developed independently using different approaches (3D finite difference, distributed model at Braunton Burrows and 1-D recharge distributed spatially at Ainsdale) to demonstrate that modelling provides a better understanding of the systems and assists in identifying what data are crucial to produce functioning groundwater models. In the paper we will assess the ability of groundwater models to represent dune system hydrology validated by historical data in order to evaluate if such modelling is reliable enough to have a role in influencing dune management.

Case studies

Braunton Burrows

Background

Braunton Burrows is a spit dune system in the south west of Britain (Fig. 2a), about 5 km long (north south) and 1½ km wide. The area comprises a series of north to south oriented dunes and slacks shaped by the prevailing on-shore winds. The slacks are mainly dry with ephemeral wet weather pools and these are sensitive to the groundwater elevation (Davy et al. 2010). The dunes system rests on an estuarine clay layer which forms the base of a small rain-fed sand aquifer (Burden 1998). The clay is underlain by Pleistocene beach deposits and/or by the Devonian Pilton Mudstone Formation which forms a low-yielding aquifer, with little to no hydraulic contact between the Devonian bedrock and the sand aquifer.

Fig. 2
figure 2

Conceptualisation of (a) groundwater flow in the Braunton Burrows aquifer and (b) model boundaries of the numerical model (displayed as 1 km grid)

Full size image

Stratford et al. (2013) have found that water levels in the dunes have been gradually declining since 1966 by about 0.7 m in winter and slightly more in the summer and this has raised concerns of the ‘drying out’ of the dune system (English Nature 1992). A groundwater modelling study was instigated to better understand the controls on the dune groundwater system and to investigate possible reasons for the observed water level decline.

Conceptual model

Water level time series data have been collected at various locations on Braunton Burrows since 1966 (Fig. 3, open cirles). Over time, some of these dipwells fell in to disrepair, and additional dipwell transects were installed in 1972 (Fig. 3, black circles) and 1992 (Fig. 3, black triangles) and have been monitored on a monthly basis ever since. This data provided the basis for the conceptual model of the dune groundwater system in Fig. 2a (Robins 2007). The dune system is predominantly rain-fed, with possible inflows from the Saunton Downs in the north. The water table is mounded within the sands with a broad groundwater divide to the east of the dunes system’s central axes (N-S), suggesting that the Burrows drain largely towards the coast. A small area of Braunton Burrows drains to the east towards Braunton Marsh where the discharge is intercepted by the West Boundary Drain with flows typically <0.5 l s−1 during the summer months (Robins 2007). The long term (1961 to 2012) average annual net (rainfall – actual evaporation) is 453 mm (as calculated from MORECS data). This suggests that in most years, rainfall recharge sustains the groundwater system.

Fig. 3
figure 3

Location of current (full circle, triangle) and past (open circle) water level monitoring wells at Braunton Burrows

Full size image

Model development

Model development consisted of a series of steps including:

  1. (1)

    Grid construction and initial parameterisation

  2. (2)

    Sensitivity analysis

  3. (3)

    Dynamic balance simulation

  4. (4)

    Model calibration.

Grid construction and initial parameterisation

A simple one-layer groundwater flow model was constructed in ZOOMQ3D (Jackson and Spink, 2004) using a uniform grid (cell size 200 × 200 m) (Fig. 2b). Boundary conditions were assigned based on the conceptual understanding provided by Robins (2007). The western and southern model boundaries were located half-way between the low and high (tidal) water mark and constraint by a Dirichlet-type head boundary, representing the mean sea water level. In the north and north-east, no-flow boundaries were assigned representing the contact to the Devonian cliffs and marshland clays, respectively. The south-eastern model boundary coincided with the Western Boundary Drain (WBD) and was represented by series of one-way leakage nodes (Cauchy-type). The model is spatially-distributed, predicting the spatial distribution of water levels within the dune system as a function of time and a variety of spatial variables (i.e. recharge, aquifer depth, permeability).

Previous studies suggest that there is little spatial variation in aquifer depth, particle size or hydraulic conductivity across the dune system (Burden, 1998; Allen et al., 2014). Hence, a uniform distribution of hydraulic conductivities (10 m d−1) was initially assumed across the entire model domain. The aquifer base was set to 0 m AOD, based on evidence presented by Burden (1998) who identified the presence of an impermeable, silty/clay layer at approximately that depth (mean sea level) at five locations across the site. Recharge was calculated using the recharge model ZOODRM (Mansour and Hughes 2004) with rainfall and Evapotranspiration data from the Bideford weather station (10 km inland) and from MORECS (40 km × 40 km integrated value), respectively. The model was run from 1st January 1970 to the 31st December 2014 using a daily time step. The calculation method used is the modified FAO (Hulme et al., 2001) using soil parameters from soil maps and crop parameters from a combination of Land Cover Map 2000 (Fuller et al., 2002) for their spatial distribution and literature values (Hulme et al., 2001). Temporally varying recharge was calculated, but the same value was applied across the entire model domain, i.e. assuming that there is no spatial variation in recharge across the site.

Sensitivity analysis

To examine the sensitivity of the model results to the input parameters and model geometry, a series of steady state model simulations were performed. A base run was carried out by running the model using long term average (LTA) recharge inputs until it reached steady state (after about 5 years). Model sensitivity was then tested by modifying a single model parameter and assessing the resulting change by comparing the outputs to those from the original base model as well as to long term average groundwater levels observed at the 17 active monitoring wells across the dune site (see black circles and triangles in Fig. 3). Sensitivities to hydraulic conductivity, aquifer base level, boundary location, leakage node elevation and leakage coefficients were tested as detailed in Table 3.

Boundary conditions were tested and revised as part of the sensitivity study. To address uncertainty around the north-east boundary, a field visit was conducted to map relevant drainage features. The survey revealed the presence of a small stream, entering the model area from the north, as well as a small drainage ditch along the north-eastern boundary. Subsequently, features were included in the model by adding river nodes as well as one-way head dependent leakage nodes along the northeastern boundary (Fig. 4a).

Fig. 4
figure 4

Sensitivity analysis: (a) refined model grid (200 m x200m) and boundary conditions (small red circles = no flow, small green circles = constant head, small blue circles = river, big green circles = one way leakage), (b) simulated groundwater head distribution (uniform parameter distribution)

Full size image

From the sensitivity analysis, a best fit model was selected and adjusted within the tested parameter range (Table 1) until the modelled groundwater heads (e.g., Fig. 4b) agreed reasonably well with observed water levels. Agreement was tested by comparing modelled (steady state) water levels against observed LTA water levels for three transects across the site (see Fig. 3). The model which produced profiles that best matched the envelope of minimum-maximum observed water levels was selected as the best fit model for further simulations (Fig. 5).

Table 1 Parameter ranges and model geometries tested as part of the sensitivity analysis
Full size table
Fig. 5
figure 5

Identification of best fit model through comparison of modelled LTA water levels along transect 1, 2 and 3 (shown in Fig. 3) against minimum and maximum observed water levels along the transects

Full size image

Dynamic balance simulation

A dynamic balance simulation was performed to obtain the correct initial conditions for the time-variant simulations (Rushton and Wedderburn 1973) and to ensure that the model reproduced the seasonal, time-variant dynamics of the dune system adequately. It consisted of a series of repetitive model runs applying monthly average recharge inputs over a cycle of three years. Successive runs were identical except for initial heads, which were replaced by the computed heads from the preceding run. Outputs were compared to the previous run in the series as well as to observed LTA monthly water levels. The runs were repeated (each time with a new initial head condition) until the computed heads followed the same cyclical pattern during consecutive runs (i.e. heads from the latest run (blue line in Fig. 6) were identical to those of the previous run (red line in Fig. 6)). A dynamic balance was reached when the net change in computed heads over a year was zero, despite changes in aquifer storage between the months. The dynamic balance head distribution provided the initial groundwater heads for all subsequent simulations.

Fig. 6
figure 6

Final dynamic balance run for selected observation well (red = previous model run, blue = latest model run) and observed LTA water levels (green)

Full size image

Model calibration

Model calibration was carried out using time-variant simulations. Since monthly time steps could not represent the observed system dynamics accurately, a daily time step was used. The ZOODRM model was re-run (as detailed previously) employing daily rainfall data from the CERF model - a regionalized rainfall-runoff model (Griffiths et al. 2008) which includes a national data base of daily rainfall (measured by the UK Meteorological Office) at the 1 km scale. Rainfall data from CERF were available for the period from 1983 to 2006, determining the modelling period for all subsequent runs.

Calculated head distribution from the dynamic balance was used as initial head conditions in all calibration runs. The model was run at a daily time step for the period from 1983 to 2006. Outputs were compared against observed water level data for the 17 observation wells (black circles and triangles in Fig. 3), visually by plotting modelled against observed water levels and quantitatively, by calculating the Nash-Sutcliffe efficiency (NSE) and the Root Mean Square Error (RMSE) for each run.

Calibration started with a simple base case of uniform transmissivity distribution (i.e., uniform hydraulic conductivity (K) and aquifer base elevation (z)) and porosity across the entire model domain. Model complexity was gradually increased during the calibration process to include 4, 6 and 12 zones across which transmissivity was varied by changing K or z. Initially, delineation of zones (4 & 6 zones) was based on the understanding that (1) the bedrock surface beneath the dunes, and hence the base of the aquifer (z), is sloping towards the sea (May 2001) and (2) that the hydraulic conductivities (K) change inland towards Braunton Marsh (Burden, 1998). Later, a pilot geophysical survey was initiated using passive seismic (TROMINO®) technology to better constrain aquifer base elevations across the model area. Using these data, a more detailed delineation of transmissivity zones (12 zones) within the dune system was derived.

The use of transmissivity zones was preferred in this study to reduce the number of ‘free’ parameters for calibration (Refsgaard 1997). Instead of varying aquifer base elevation (z) and hydraulic conductivity (K), one parameter, usually z, was fixed for the entire study area while the other one (K) was adjusted during calibration to simulate variations in transmissivity (T). Hence, values of K derived in this study present an adjustment factor for T, not a hydraulic conductivity per se.

Monte Carlo runs were conducted for both, simple (1 zone) and more complex (4, 12 zones) transmissivity distributions to identify the best parameter sets for calibration. Each model setup was run for up to 10,000 iterations, varying transmissivity (K or z) or aquifer storage (S) parameters during each iteration and for each zone randomly within a predefined range (Table 2). A variety of error measures (including NSE) were computed. However, from the computed error measures (which included NSE) it was not possible to identify best parameter values for the different zones (e.g., Fig. 7). Hence, a trial-and error method (Refsgaard and Storm 1996) was adopted, manually adjusting either K or z or S within the different zones (while keeping all other parameters constant) and observing the resulting model fit visually through plots and quantitatively through calculation of the NSE and RMSE.

Table 2 Parameter ranges for the Monte Carlo simulations
Full size table
Fig. 7
figure 7

Monte Carlo simulations: (a) Transmissivity zone delineation used for MC runs, (b) Impact of varying hydraulic conductivity (K) in zones 1–4 on model fit (Nash Sutcliffe Efficiency) in OBH 1. Red circles indicate NSE values >0 (i.e acceptable model fit)

Full size image

Model validation, as commonly used in surface water modelling, is not undertaken here as it is not considered useful or even possible within the context of groundwater modelling (Anderson et al. 2015, Konikow and Bredehoeft, 1992). Because of the large numbers of parameters involved, Doherty and Hunt (2010) advise using all available data in the calibration process as omitting any data is likely to result in a poorer calibration.

Results and discussion

Figures 5- 11 and Table 3 summarise key outputs and model performance measures from the different stages of the model development process.

Table 3 Summary of key model runs: zone delineation, parameterisation and model performance measures (NSE, RMSE)
Full size table

Sensitivity analysis: The final run of the sensitivity simulations, i.e. the selected best fit model, is shown in Figs. 5a and b. Figure 5a shows a satisfactory match between the modelled water levels (solid line) and the observed minimum and maximum water levels (broken lines) along the three transects (shown in Fig. 3), with better agreement in the northern part of the dune system (Transect 1) compared to the south. Water level contours produced by the model (Fig. 5b) were very similar to observed water level contours, e.g., as presented by Burden (1998) and Willis et al. (1959b) and hence, this model was taken to be a satisfactory starting point for further model development.

Dynamic balance: Results from the dynamic balance simulation are shown in Fig. 4 for selected observation wells. Initial runs (not shown) displayed a time lag of around one month between the overserved and the modelled seasonal signal, representing the average time needed for effective precipitation at the surface to reach the groundwater table (Moon et al. 2004). Applying a lag correction of one month to the recharge data resulted in better agreement between modelled and observed system dynamics (Fig. 4).

Model calibration: Fig. 8 shows output from the Monte Carlo simulations for the uniform (1 zone, see Table 3) transmissivity distribution. In this example, the base of the aquifer was fixed at −2 m AOD, but transmissivity of the zone was varied by randomly varying hydraulic conductivities within the range defined in Table 2. In Fig. 8, K values for each run are plotted against the corresponding NSE (green dots representing NSE ≤ 0 and red dots representing NSE > 0) for the 17 observation wells (OBHs) in the study area. The plots show that the K values required to achieve a good model fit vary between the different monitoring wells, suggest that transmissivity (T) varies across the model area. There is a general decrease in K (and hence T) from the western to the eastern boundary, which could indicate a decrease in aquifer saturated thickness (e.g. due to a sloping bedrock surface (May 2001)), a decrease in hydraulic conductivity (Burden,1998) or both. At a number of observation wells, including OBH 6, 12 and 17, NSE remains below zero, i.e. water level fluctuations at these wells could not be simulated satisfactorily by any of the runs. Their location near the edges of the model suggest that the boundary conditions may need to be revised in order to improve the model fit in these areas. A model run was conducted using a K of 10 m d−1. Model parameterisation and performance (NSE and RMSE) are recorded in Table 3 for the different wells. Model fit is predictably poor with an average RMSE of 1.44 (range 0.34–3.29) and a NSE > 0 at only 2 (12%) of the 17 observation wells. This confirms the conclusion that a single T zone (z = −2 m AOD, K = 10 m d−1) cannot describe the dune system hydrology adequately.

Fig. 8
figure 8

Monte Carlo simulations: (a) Location of observations well and delineation of transmissivity zone, (b) Plots of hydraulic conductivity (K) against Nash-Sutcliffe efficiency (NSE) for all observations wells (arranged along transects 1–3)

Full size image

An attempt was made to use Monte Carlo for an increasing level of complexity, e.g., including 4 transmissivity zones (Fig. 7a) within which K was simultaneously varied within a predefined range (Table 2). The results are shown in Fig. 7b for OBH 1, and demonstrate that this introduced too many degrees of freedom and resulted in non-identifiability. Consequently, calibration was continued by manual parameter assessment through a number of simulation runs.

Defining four and six transmissivity zones (as shown in Table 3) to simulate a gradual increase in transmissivity towards the sea improved overall model performance. Both models were able to reproduce the general water level distribution (i.e. mounding of the water table to the east) Fig. 9 as well as the seasonal dynamics of the dune groundwater system Fig. 10. Model outputs from the 6-zone model compared reasonably well with observed water levels in the central part of the dune system (e.g., OBHs 1–5, 11, 15) (Fig. 10, Table 3), but water levels were overestimated (e.g, OBHs 6, 8, 12, 13) or underestimated elsewhere (e.g., OBHs 9, 10). There is a noticeable discrepancy between modelled and observed water levels in all the observation wells. This discrepancy is greatest during the last three years of the simulation (2003–2006) (e.g., OBH 11, Fig. 10). It cannot be explained by model performance alone but is more likely to be linked to the driving data (i.e. rainfall) as will be discussed later. Consequently, the error measures in Table 3 were calculated for the period 1983–2003 (i.e., excluding the last three years of simulation) to account for this.

Fig. 9
figure 9

Water level contours generated by (a) the 4-zone and (b) the 6-zone model

Full size image
Fig. 10
figure 10

Model output from the 6-zone model: (a) Location of observations well and delineation of transmissivity zones, (b) Plots of observed (black circles) versus modelled (green line) water levels for all observations wells (arranged along transects 1–3)

Full size image

Overall model performance (Table 3) was better for the 6-zone model than for the 4-zone model, increasing the number of observation wells with NSE > 0 from three (18%) to five (30%) out of a total of 17. The RMSE also improved from an average of 0.94 (range 0.36–2.36) to 0.64 (range 0.32–1.36). Through adjustment of K, model fit could be improved locally but this caused deterioration in model fit elsewhere, hence the overall model fit remained poor. This indicates that the spatial distribution of transmissivity (and probably also recharge) is more complex than represented by these models. In the absence of additional data, it was felt that calibration could not be improved without making unjustifiable assumptions and that further data, relating to aquifer base topography or K distribution were required.

A pilot geophysical survey of the bedrock depths was initiated as part of this modelling study. Based on the data derived from this work, a more detailed delineation of transmissivity zones was possible. The findings that are of relevance here, include the following points: (1) the bedrock topography underlying the model area is more complex than initially assumed, (2) the depths to bedrock generally increases from around -10 m AOD in the west to 0.5 m in the east (confirming that the underlying Devonian shelf platform is tilting seawards as suggested by May (2001)) with a ridge of slightly higher relief in the central part of the dune system, (3) Depths to bedrock of −34 m to −17 m AOD are observed at the northwestern boundary of the dune system and at locations in the northeast and southeast of the dune system, respectively. The extent of these features needs to be investigated further, but they may represent deep channels or depressions D. Boon. 2015 (personal communication), similar to the Pleistocene ‘rock-valleys’ near Plymouth (Eddies and Reynolds 1988) or the ‘buried channel’ in Bideford Bay, immediately south of the Braunton Burrows site, which was estimated to reach a depth of about −30 m AOD (McFarlane 1955).

Translating the new understanding into a transmissivity distribution resulted in the definition of 12 transmissivity zones (Table 3, Fig. 11), including a high transmissivity zone in the north (possibly representing a sediment-filled channel). The bedrock depressions along the eastern boundary were not considered as separate zones in this distribution as the shape and extent of these features could not be sufficiently defined from the available data. A new aquifer base of −5 m AOD was assigned, representing the average bedrock depths observed during the geophysical survey. Monte Carlo runs were conducted for each zones separately to provide initial K estimates, which were then adjusted using a trial-and-error calibration method.

Fig. 11
figure 11

Model output from 12-zone model: (a) Location of observations well and delineation of transmissivity zones, (b) Plots of observed (black circles) versus modelled (green line) water levels for all observations wells (arranged along transects 1–3)

Full size image

The results are summarised in Table 3. It shows a considerable improvement in overall model performance, with 12 from 17 observation wells (70%) displaying a NSE > 0 compared to 5 from 17 (30%) for the 6 zone model. Improvement in model fit is also reflected by the lower average RMSE of 0.40 (range 0.25–0.84) compared to the previous 0.64 (0.32–1.36).

The water level dynamics are generically well reproduced by the model (Fig. 11), except in OBH 6, where water level variations are less dynamic than predicted. Overall water levels are predicted correctly by the model, except at OBH 10 (not shown), where water levels are underestimated by about 50 cm, and OBH 17, where they are overestimated to a similar degree.

There is still a discrepancy between modelled and measured data, which is observable in all hydrographs in Figs. 10 and 11. Examining these hydrographs more closely shows that this discrepancy appears to start in the early 2000s. For example, the observed groundwater levels show a significant fall in late 2001 before levelling off, whereas the modelled response is a continued cyclical decrease in groundwater head. There appear to be a number of rainfall events (smaller and one larger) in the observed data which are not reproduced in the modelled response. For example, the observed response shows an event in early-mid 2004 which leads to a rise in groundwater heads in all OBHs across the site but is not seen in the modelled response. This event could result in the groundwater system receiving a large mass of recharge which slows the year on year reduction of groundwater heads within the dune aquifer. If these events were not captured in the gridded rainfall data, used here to drive the model, this could explain why modelled water levels decline whereas the observed data show a level response.

Ainsdale (Sefton coast)

The Sefton Coast stretches 20 km between Crosby and Southport in Merseyside, UK. The coastline consists of a dune system 2–3 km wide, running parallel with the coast (Fig. 12). Average rainfall in the area is approximately 800 mm y−1 and evapotranspiration is around 500 mmy−1. The dune system comprises approximately 25 km2 of open dunes, forest and golf courses and the built up areas of Formby, Ainsdale and Southport. The dunes rise to 20 m above sea level and rest on a bed of peaty clay which effectively seals it from deeper Triassic sandstone aquifers. The dune system is home to a number of nature conservation sites and Nature Reserves, the largest of which is Ainsdale Sand Dunes National Nature Reserve, established in the 1960s. The Reserve aims to protect ecosystems which include the seasonally flooded low-lying areas (slacks) between the dune ridges.

Fig. 12
figure 12

Ainsdale Nature Reserve on the Sefton coastline. Arrows indicate directions of groundwater movement from the high point of the groundwater mound

Full size image

Conceptual model

A set of 12 observation wells were installed at 500 m spacing on the dune system in the early 1970s to monitor groundwater levels on a monthly basis. Contour maps of the levels demonstrated the existence of an elongated mound of groundwater which runs parallel to the coast and from where groundwater drains inland and towards the coast as shown in Fig. 12. The highest point of the groundwater mound is 10 m AOD and the typical seasonal variation in groundwater levels is +/− 0.5 m (Fig. 13). Highest water table levels occur at the end of the winter season around March when many of the dune slacks are flooded to a depth of 0.1–0.3 m (Clarke and Sanitwong Na Ayuttaya 2010). There appear to be sequences of wet years and dry years in the data but no long term trend can be detected from this 40 year record (Fig. 13).

Fig. 13
figure 13

Monthly observations of water table levels at Ainsdale 1972–2013. Red line = ground level. Shaded bands are based on Ranwell (1953) classifications of slack type – wet/transitional/dry

Full size image

Model development

Clarke (1980) and Clarke and Sanitwong Na Ayuttaya (2010) describe groundwater models based on a water balance between rainfall, evapotranspiration and groundwater flow (Fig. 14). This represents a point/lumped model to predict water levels at a selected location within the dune system as a function of time (but not of space). Key inputs are rainfall and an estimate of actual evapotranspiration. Rainfall data were available from a local meteorological site, which also provided the climatic data to estimate potential evapotranspiration using the Penman Montieth equation (Allen et al. 1994). Actual evapotranspiration is calculated using a soil root zone storage model and assumed that surface runoff occurs rarely and only when multiple interconnected slack floors are flooded. When the soil water in the root zone is full (i.e. the soil moisture deficit is zero), and any surplus water is assumed to percolate vertically down to the water table. The groundwater balance is computed on a daily time step. Groundwater flow out of the dune system is calculated using the Dupuit formula. The flow is a function of the sand permeability (K) and thickness (b) and the water table gradient towards the sea or inland. In this model the permeability and thickness were assumed to be constant over the study area. At the end of each day, the elevation of the water table at the modelled location is updated based on the balance between the groundwater outflow and the aquifer recharge. On dry days the soil is assumed to dry out first and the water table will descend based on its current height and the dune permeability. On very wet days surplus water percolates to the water table and raises it by a vertical distance of recharge (in mm)/n where n is the effective porosity of the dune sand. When the water table is above ground level, n is set to 1. The model then carries the computed water table level forward to the next day (time step) when the calculation is repeated. This produces a time series of modelled water table levels.

Fig. 14
figure 14

Conceptual model of the components of the water balance in the dune system. The water table (red line) slopes towards the sea and intersects low lying areas to create humid or flooded slacks

Full size image

Sensitivity analysis and model validation

The model was initially developed for one site, located in the open dunes in northern part of the area. It was calibrated using parameter optimization based on measured values of sand permeability (K) and dune system thickness (b), the spatial distribution of vegetation/open water areas (which affect evaporation rates), and soil water storage characteristics such as dune total porosity (n) and the effective porosity (ne), which changes from 0.28 in clean deeper dune sand to 0.2 in the vegetation close to the surface and 1.0 above the surface. The model was run for a period of 30 years and with successive refinements to values for the dune sand permeability, porosity and aquifer thickness. Typical values for K for blown sand are reported between 8 m d−1 to 15 m d−1 (Freeze and Cherry 1979) and n between 0.35 and 0.4. Model results for the two extremes of K are shown in Fig. 15a and b. Setting the permeability to the upper value of 15 m d−1 results in too much leaving the model as groundwater flow, and the predicted water table levels are consistently lower than the observed levels (Fig. 15a). When K is set to the lower value of 8 m d−1, not enough water leaves the system as groundwater flow and the predicted levels are higher than observed (Fig. 15b). Over the 30 year test period, it was found that a value of K = 11 m d−1 resulted in the average modelled water levels being closest to the observed data (Fig. 15c). The amplitude of seasonal variation in modelled water table levels was found to be highly dependent on the effective porosity. The laboratory determined dry dune sand porosity was 0.39 but model outputs suggest the amplitude between winter and summer was too small. As the natural dune system retains some water in capillary spaces between the sand grains, the effective porosity (ne) of the sand should be <0.39. A sensitivity analysis showed that when the water table was >0.5 m below ground level the best fit value of effective porosity (ne) was 0.28. This was adjusted as the water table rose closer to the surface and in the top 50 cm ne was set to 0.2 and above the ground level, ne = 1.0 (i.e. 100% air space).

Fig. 15
figure 15

Model validation for water table levels at Ainsdale. Sand dune properties, K (hydraulic conductivity) and n (effective porosity), were successively adjusted until the best fit was obtained between the model and the observations. Note that the aquifer thickness (‘b’) was kept constant at an assumed thickness of 60mbsed on BGS local borehole logs

Full size image

Two measures of goodness of fit were used to evaluate model performance – the root mean square error between observed and modelled levels (RMSE) and the Nash-Sutcliffe efficiency coefficient (NSE). Values of these measures are shown in Fig. 15 when the model is applied to a typical observation well in the dune system. The best model fit (RMSE = 0.12 m, NSE = 0.84) was found to occur when permeability (K) was set to 11 m d−1, effective porosity (ne) to 0.28 for the dune sand deeper than 50 cm, rising to 1.0 at the surface (i.e. open air) and the dune sand thickness (b) equal to 60 m (Fig. 15d).

Following testing and calibration, the model was able to reproduce the vertical changes in water table levels within +/− 10 cm over a 30 year period. The model was applied to two other wells, one in the open dune fields close to the sea and another under the pine forest further inland. In both cases, the model parameters required adjustment and re-calibration as the vegetation cover differed. In the case of the dune fields closer to the sea, there is a larger area of slacks and fewer dune ridges. This affected the soil water balance and in winter time, the flooded slacks joined up and were able to generate some overland flow in addition to the groundwater flow. Under the pine trees, the effects of rainfall interception had to be included in the model, based on the algorithms suggested by Calder (2005). In both cases, the adjusted models were able to replicate observed water table levels. However, the general applicability of the water balance approach to different parts of the dune system required fine tuning to the specific sets of conditions at each location.

Discussion

In the previous sections, we have described two independent approaches for modelling dune groundwater systems, (1) a three-dimensional, numerical groundwater flow model (ZOOMQ3D) with spatially distributed parameters for dune system properties and (2) a one-dimensional single-point / lumped parameter model, to simulate groundwater flow in small sand dunes systems at Braunton Burrows and Ainsdale, respectively.

Developing the groundwater flow of the Braunton Burrows dune system has shown that a higher degree of model complexity/heterogeneity has been required than was initially anticipated based on the conceptual model. Whilst the system itself is of relatively modest size (8 km2) and has been subject to a variety of surveys and studies (e.g. Mcfarlene 1955; Kidson & Carr 1960; Kidson et al. 1989; Burden 1998; May 2001; Stratford et al. 2013; Allen et al. 2014), quite a high degree of uncertainty was discovered during model development relating to the definition of the model boundary conditions (e.g., drainage conditions in the northeast of the model domain) and the transmissivity distribution within the model domain.

Simplifying assumptions had to be made during the early stages of model development, but these were gradually adjusted using sensitivity analysis. The importance of such an analysis in defining the correct model boundary conditions has been demonstrated by Hunt et al., (1998).

Targeted field investigations were an important part of the model development process, specifically the geophysical survey, which provided an improved characterization of the bedrock topography and led to a more detailed delineation of transmissivity zones within the model domain. Constructing models with increasingly more complex spatial distributions of transmissivity within the model domain, led to considerable improvements in the fit between predicted and observed heads.

Calibration of transmissivity (T) (= one parameter) rather than hydraulic conductivity (K) and depths to aquifer base (z) (= two parameters) was chosen here to reduce the degrees of freedom or “free parameters” in model calibration, as proposed by Refsgaard (1997), and to reduce the problem of equifinality or non-uniqueness in numerical modelling (Beven 2006), which says that there are many different conceptualizations of a numerical groundwater flow model that fit the observed data equally well. This problem has been demonstrated by the Monte Carlo analysis (Fig. 7), where many different combinations of T resulted in the same or similar NSE values (despite the reduction in free parameters). This non-uniqueness is also the reason for why numerical groundwater models cannot be truly validated (Konikow and Bredehoeft, 1992) and why their use must be constraint by their inherent uncertainty (Wondzell et al., 2009).

In the Braunton Burrows model, model uncertainty needs to be reduced and the problem of equifinality must be further investigated. While the topography of the aquifer base (i.e. z) is relatively well-constrained by the data from the geophysical survey, there remains large uncertainty relating to the distribution of hydraulic conductivity (K) within the dune system. The inability of the model to reproduce the observed water level response in some parts of the model may indicate that the distribution of hydraulic conductivities across the model domain is not as uniform as Burden (1998) and Allen et al. (2014) suggest. Recent drilling in the central part of the dune system has revealed the presence of coarse sands and gravels at depths within the dunes (D. White, 2015, personal communication), hence challenging the general assumption that sand-blown dune systems consist of uniform, well-sorted granular sediments throughout. Furthermore, the lack of water level data in the north of the dune system must be addressed and an improved understanding of how rainfall events in 2004 influence recharge is required.

Groundwater modelling is often considered to be a linear process, consisting of a series of successive steps as, for example, outlined in Hobma et al. (1995). The present study illustrates that it should be a more iterative process and that modelling should be carried out in conjunction with field studies, i.e. not as a postscript, but as ongoing interaction. This interaction is required throughout the investigation process. The importance of feedback between numerical modelling, conceptual modelling and data collection is increasingly recognised in the modelling of dune systems (e.g. Spring (2005)) as well other groundwater systems (Voss 2011).

The simple, conceptual modelling undertaken for the Ainsdale Dune system demonstrated that such a model could be calibrated with reasonable generic hydraulic parameters (K = 11 m d−1; ne = 0.28). There are a number of observation boreholes on the Ainsdale site and the model transferability was tested by applying it to two more boreholes. However, adjustments were required as at one location the dune slack generation was more pronounced, and the other was covered by pine trees. For the latter, the recharge module of the model needed adjustment to take into account interception. The changes required by the model illustrate the variability within a small site. The subtle variation of land use and surface response required changes to the model structure thus reinforcing the need for interaction and subsequent flexibility within the modelling methodology.

The modelling process is often represented by a flow chart to illustrate the interaction between different activities. Figure 16 shows an extension to those typically presented in the literature (e.g. Anderson et al. 2015) and illustrates the approach adopted in this study. It shows that modelling is not just a passive activity, but one that involves developing a conceptual understanding of the system, encapsulating this into a numerical model and then testing the model results against observed data (e.g., during model calibration). This process enables the observed data to be actively interrogated and not just examined passively. The result is that the model can show the user whether their conceptual model, when translated into a groundwater model, reproduces the observed response.

Fig. 16
figure 16

Representation of the groundwater modelling process

Full size image

In first runs of models, a poor fit is usually found between the water table levels observed in the field and the model-generated data. This often forces the modeller to return to the conceptual model or to question the data used to parameterise or run the model (e.g. the representativeness of the available dune property data or the scale at which rainfall or evapotranspiration data are employed). This interactivity is key to heuristic learning, i.e. learning by someone doing something themselves, and helps guide the user through developing their understanding of the system and identifying important controls. Voss (2011) points outs that rather than use automatic calibration the interaction between the user and the model should be “learning something from a model-based analysis”.

Conclusions and recommendations

  • Data collection in dune systems was initially instigated in the UK in the 1970’s to monitor system behaviour, but managers of these systems realised that large variations in water table levels could not be explained as a simple response to rainfall or changes in land use practices. Groundwater models provide effective tools to explore what controls such water level variation and to investigate the hydrology of dune systems.

  • Two modelling approaches were presented in this study which differ in complexity, but which both have produced models that successfully describe the water table variations in the simulated dunes systems.

  • A high degree of model complexity and between-site adaptations in model structure were required to be able to reproduce the dune systems’ behaviour. This illustrates that these systems are hydrologically complex despite their small size and clear boundaries. Factors such as bedrock topography, land use and local climatic conditions are important controls on their hydrological response and hence, need to be described adequately by the model. The general assumption that sand-blown dune systems have a uniform and well-sorted sediment composition, and hence a uniform K distribution, needs to be investigated for each dune site individually and treated cautiously until proven by representative field surveys.

  • Because of their small size, models representing such small aquifer systems are more sensitive to uncertainties in input data (i.e., basic driving data such as rainfall and evaporation), model geometry (positioning of model boundaries) and model parameterisation as well as to changes (i.e. land use) that occur during the simulation time period. Thus, it is critically important to collect time series data (e.g., rainfall, evapotranspiration, water levels) at relevant temporal and spatial scales and also to explore the dune system for variations in thickness, land cover, grain size distribution and permeability.

  • The case studies have demonstrated that modelling can support the development of systems understanding through repeated interrogation of the model outputs, observed data and conceptual understanding. This process usually raises question relating to conceptual understanding and data availability. It is therefore important that modelling is carried out in parallel to field investigations and monitoring, ideally during earlier stages of the research project, so that it can inform researchers and managers of important data needs or identify locations where additional monitoring or data collection would result in improved system understanding and reduction in uncertainty.

  • Modelling should be an iterative process which guides the user towards a better understanding of the modelled system. Deficiencies in the model in reproducing the observed response, both spatially and temporarily, point the modeller to where further data and/or process understanding is needed.

  • Given that the role of modelling is to raise questions as well as to answer them, this study has demonstrated that this applies even in small systems that are thought to be well understood.