Interactive comment on “On the validity of modeling concepts for (the simulation of) groundwater flow in lowland peat areas – case study at the Zegveld experimental field” by P. Trambauer et al

The work by Trambauer et al. if focused on the possibility of reliably simulating the flow field in the Holocene peat layer that outcrops in a large portion of The Netherlands. This task is accomplished by both experimental and modeling activities. The objective is clearly very interesting from the scientific point of view as the hydrologic response of peats is well known to be very complex. The goal is also very important from an applicative point of view because this peat unit constitutes the seal for the underlying sandy aquifer.

leads to transient behavior.However, using transient formulae, the slug tests conducted for different initial groundwater heads showed that there was hardly any evidence of a variation of the hydraulic conductivity with the applied head differences.Therefore, Darcy's law can be used for typical peat soils present in The Netherlands.
The heterogeneity of the subsoil and the apparent applicability of Darcy's law were taken into account for the detailed heterogeneous model that was prepared for the research area.A MODFLOW model consisting of 13 layers in which 4 layers represent the heterogeneous CCL was set up for an average year, assuming steady-state conditions; and for the winter of 2009 to 2010, adopting transient conditions.The transient model was extended to simulate for longer periods with the objective of visualizing the flow paths through the CCL.The results from these models were compared with a 10 layer model, whereby the CCL is represented by a single layer assuming homogeneity.From the comparison of the two model types, the conclusion could be drawn that a single layer schematization of the CCL produces flowpath patterns which are not the same but still quite similar to a 4 layer representation of the CCL.However, the single layer schematization results in a considerable underestimation of the flow velocity, and subsequently a longer travel time, through the CCL.Therefore, a single layer model of the CCL seems quite appropriate to represent the general flow behavior of the shallow groundwater system, but would be inappropriate for transport modeling through the CCL.

Introduction
In the western part of The Netherlands there is a great need to understand the flow of groundwater in lowland peat areas and its interactions with surface water.Generally in these areas an aquitard consisting of peat and clay overlays a sandy Pleistocene aquifer.The aquitard or semi-permeable layer is often referred to as a Complex Confining Layer (CCL) due to its lithological and hydraulic heterogeneity (Dufour, 2000;Bierkens, 1994;Weerts, 1996).This layer, which covers a large area of The Netherlands, plays an important role in the interactions of surface water and groundwater systems in the low lying parts of the country, such as protecting groundwater from pollution.
The groundwater flow models currently used in peaty lowlands are likely a too simplified representation of the hydrological reality because of the following two reasons: 1.The CCL is often represented as a one model layer whereby different parameter values in this layer are lumped together.In another approach, the CCL may be considered as a single peaty model layer representing a phreatic aquifer resting on top of a Pleistocene sandy aquifer with horizontal flow.The clayey components in the peat representing aquitards are modeled as a model resistance function between the peaty aquifer and the underlying sandy aquifer.However, for the CCL in particular, it is thought that a more complex representation is needed whereby a multi-layer set up takes better into account the heterogeneity of the clayey peat layer.
2. In many studies, groundwater model codes such as MODFLOW (McDonald and Harbaugh, 1998) are used which are based on physical concepts that use a combination of Darcy's law and the mass balance equation for the computation of groundwater heads and flows.However, according to Ingram et al. (1974) and Waine et al. (1985), the flow in humidified peat does not follow Darcy's law.The questionable applicability of Darcy's law in peat areas has an immediate effect on the modeling tool that can be used for management tasks.Studies suggesting the non-validity of Darcy's law for peat soils indicate that the flow differs from the flow in a rigid mineral soil.An increase of the hydraulic conductivity of the peat when the groundwater head gradient increases, has been observed (Ingram et al., 1974;Waine et al., 1985).However, the applicability of Darcy's law and the robustness of the hydraulic conductivity for peat soils apparently depend on the degree of decomposition or humification of the peat.
The main objective of this paper is to investigate whether the described reasons lead to an inadequate modeling of groundwater flows in the CCL.Geo-hydrological data were collected and analyzed at an experimental field within a research area located in a typical peaty lowland area in the vicinity of Zegveld, The Netherlands.Data were collected to investigate the heterogeneity of the peaty CCL and to verify the applicability of Darcy's law in peat soils.Subsequently, groundwater models with a complex heterogeneous or simplified homogeneous CCL were built.The models were compared with each other to assess the effects of simplifying the CCL and to test the applicability of Darcy's law in peat.

Research area
The research area was located in the Zegvelderbroek polder at approximately 2.5 km north of Zegveld, a small town in the province of Utrecht.Within the research area, experimental fields to carry out data collection programs have been defined.To collect the data for the current research, a field was selected having the form of a rectangular plot of land and completely surrounded by ditches (Fig. 1).The water levels in the ditches were kept at −2.95 m during the summer months and at −3.0 m during the winter months.The plot had a surface area of 1.57 ha.The topography was very flat with land surface elevations averaging around −2.30 m.The grass covered soils received a mean precipitation of 826 mm yr −1 and the potential mean evapotranspiration was estimated at around 545 mm yr −1 .The precipitation surplus of about 280 mm yr −1 was mainly recharging the shallow groundwater.
The CCL in the area is part of the Holocene layer.This CCL has a thickness of approximately 6 to 7 m.The lithology of the CCL is mainly dominated by peaty and clayey materials which were formed and deposited in an environment with slow moving meandering rivers where the peat grew in still waters at considerable distances from the main water courses, and the clays originated from flood depositions.Clays of a marine origin are also present.Sand deposits at the margins of the area occur as well and can be related to bed deposits of the ancient water courses.The Pleistocene below the Holocene consists of medium to coarse sandy deposits of considerable thickness (around 18 m).According to Bierkens (1994) and Weerts (1996), these sandy deposits, also referred to as the Kreftenheye Formation, were deposited during the Weichselian ice age and originated from braided river systems.
The research area is characterized by an extended network of ditches and canals to drain the polder areas in the winter season and to provide high groundwater tables during the summer.Thus, there is a strong interconnection between the groundwater and surface water.The organization responsible for water management in the area is the Water Board named Hoogheemraadschap de Stichtse Rijnlanden (HDSR).

Field and laboratory work
To investigate the heterogeneity and the applicability of Darcy's law in groundwater flow modeling, activities at the Hydrol.Earth Syst.Sci., 15, 3017-3031, 2011 www.hydrol-earth-syst-sci.net/15/3017/2011/  and 2).Another two boreholes were drilled at these two sites until halfway through the CCL.Detailed sample descriptions of the collected soils were completed in the field.In addition, a geophysical survey was carried out in order to have more information on the aerial consistency of the layering detected during borehole drilling.The vertical electrical sounding (VES) method following the Schlumberger electrode layout was used.Slug tests were carried out in 8 shallow boreholes with depths up to 2.2 m to verify whether the hydraulic conductivity in the peat layers is constant or varies with the imposed groundwater heads.The slug tests were completed as follows: after the borehole was drilled and cased to avoid the borehole from collapsing and the water level (head) inside the hole equaled the groundwater table, water was instantaneously extracted from the borehole with a bailer.The groundwater heads inside the boreholes were measured with a measuring tape attached to a float at 10 s intervals from the start until the level (head) was almost stable in about 7 to 10 min.
Through laboratory tests, the relation between the variability of the hydraulic conductivity and the decomposition degree of the peat layers was established.In every borehole used for the slug tests, one or two peat samples were taken at depths ranging from 0.7 m to 2.1 m with the aim to analyze the different soil types.The decomposition degree was determined according to the von Post method, the extracted carbon analysis technique, and through the measurement of the absorbance (Price et al., 2005;Wang et al., 2004;Klavins et al., 2008).
Finally, time series of the groundwater levels were collected.The two boreholes penetrating through the CCL were equipped with piezometers having filter screens at 7 m depth in order to measure confined groundwater heads in the first Pleistocene aquifer.The other holes drilled until halfway through the CCL obtained piezometers as well.Filter screens were placed in these holes at 3.5 m depth to measure phreatic groundwater heads in the CCL.Shallower piezometers installed within the frame of earlier studies were present as well.Pressure transducers set up in all the piezometers guaranteed a continuous set of time series of groundwater heads during the fieldwork period (Fig. 3).

Results of the field investigations
The soil descriptions showed a top layer of clayey and silty peat of around 1.5 m thickness.Below this layer the CCL consisted mainly of a highly moisturized, muddy material mixed with peat.In some places a lot of wood was present (Fig. 4) but there were localized places where only peat was found.Below the mud layer, a peat layer mixed with mud and some wood was observed.Near the 5 m depth mark, a thin layer of clay was located in one of the boreholes that was below a layer consisting of a mixture of clay with peat (humic clay).At the depth of 6.3 to 6.5 m the sandy Pleistocene aquifer was found.
A cross sectional view through the CCL based on the soil descriptions obtained at the deeper and shallower boreholes indicates continuity in lithology for most layers whereby the differences in layer thickness should be noted (Fig. 5).However, not all the layers were found to be continuous as is demonstrated by the soil description of borehole 1, which had a muddy peat layer below the first layer that was not present in borehole 2. Underpinning this discontinuity is the isolated layer of peat detected in borehole 3 and the thin layer of clay that was observed in borehole 1.These layers are not present at borehole 2. Although there is a consistent lithological succession from peaty units into more clayey layers, the substantial variations in layer thickness and the discontinuity of layers in a lateral sense prove that the CCL is a heterogeneous unit.
The geophysical VES measurement completed at the research site attained an investigation depth of about 20 m which is well into the sandy Pleistocene aquifer.The interpretation of the measurement allowed a 2-layer schematization of the CCL indicating layer thicknesses of 1.2 m and 6.5 m (Fig. 6).These layers are thought to coincide with the first clayey/silty peat layer and the underlying mud and peat layer.The clay and clayey peat layers could not be differentiated due to the small thicknesses of these layers and to the fact that the latter layer contained a lot of peat.The resistivity contrast between the clay containing layers and the mud and peat layers was apparently small, which did not allow the interpretation of more layers.The third layer found in the measurement represents the sandy Pleistocene deposits.The VES interpretation shows that at this location in the middle of the research site the general lithological buildup of the CCL is similar to the sites where the boreholes allowed a detailed layer description.

Validity of Darcy's law
The slug tests carried out at the experimental field to test the applicability of Darcy's law in the shallow clay, silt and peat, and muddy peat layers were first analyzed following methods based on steady-state groundwater flow assuming a      rigid system and negligible effects on the water table.The van Beers (1963) approach and the Hvorslev method as suggested by Surridge et al. (2005) were selected.The steady state formula as formulated by Surridge can be expressed as follows: where K h (m d −1 ) indicates the horizontal hydraulic conductivity, t(d) is time, h t (m) and h 0 (m) are the groundwater heads caused by the removal of water at time t and at the start of the test,respectively.A (m 2 ) is the inside cross-sectional area of the piezometer standpipe and F (m) is the shape factor for the piezometer intake (Surridge et al., 2005).The heads are the vertical distances measured between water levels during the test in the hole and the static water level in the soil.
The formula shows that a plot of ln (h t /h 0 ) against t produces a straight line since the other parameters in the expression are assumed to be constant.
Plots of ln (h t /h 0 ) or head (level) displacement against t for the research site do not produce straight lines at all sites.A curved line is shown for the slug test carried out in Borehole 5, where the tested material consisted of peat with a Schematization of the CCL based on borehole descriptions.The average groundwater level for an average y is at about 40 cm below the ground level.high degree of decomposition (Fig. 7).Several K h values can be computed for the measured head displacements used to construct a curve.This easily leads to the conclusion that Darcy's law is not applicable in peat areas since the K h varies with the head differences in the system.
Rather than attributing the phenomena to the nonapplicability of Darcy's law, the "curved line behavior" could be attributed to the peat not being a rigid framework for groundwater flow.Peat is compressible and therefore a transient behavior of the soil can be expected during slug tests.The use of the traditional steady-state formula should be avoided (Hinsby et al., 1992;Surridge et al., 2005) and instead transient formula have to be applied.Transient formulae to interpret slug tests in unconfined aquifers are hardly described in the literature, but the KGS approach (Choi, 2008;Esling and Keller, 2009) and the Dax expression, as suggested by Hinsby et al. (1992), can be considered.
The KGS approach is based on a semi-analytical model to estimate the storativity and the hydraulic conductivity of an unconfined aquifer for the transient situation.The model incorporates the effects of partial penetration, anisotropy, and the presence of variable conductivity well skins (Esling and Keller, 2009).The KGS formula is presented in Choi et al. (2008), where the conventional line-fitting procedures as used for the Hvorslev and the Bouwer and Rice methods are adapted by introducing a modified effective radius, R e .This parameter is dependent on the compressibility of the geologic formation as well as borehole geometry (Chirlin, 1989).Because the effective radius is replaced by R e , there is no difference in evaluating the hydraulic conductivity using either the Hvorslev method or the Bouwer and Rice method.R e is obtained after determining the value of the storativity for the compressible formation and a dimensionless compressibility parameter for the given hole and tested formation geometry.In the KGS model, the storativity is estimated through curve fitting.
Dax simplified and approximated the transient Cooper method used for fully screened wells in confined aquifers for partially penetrating wells and applied this method to unconfined aquifers with delayed yield, assuming horizontal radial flow to the well.The Dax expression can be described as follows (Dax, 1987): where K h (m d −1 ) is the horizontal conductivity, t(d) is time, and H t (m) and H 0 (m) are groundwater heads, following the injection of water, at time t and at the start of the test.The r c is the inside radius of the piezometer standpipe, L is the length of the filter pack, and the D(α) is a variable also containing time t.The α includes a term for the storativity S of the soil which represents the capacity for elastic storage depending on the compressibility of the material (Hinsby, 1992).The parameter also incorporates r w (m) which is the borehole diameter, and B, describing borehole geometry.
Where the unconfined aquifer is showing a delayed yield response, the S can be obtained using the Cooper method (Cooper et al., 1967).Using steady-state and transient formulae, K h values were computed for the 8 shallow boreholes using a spreadsheet (Table 1).Head (level) displacement against time t plots were prepared for all the holes.Five holes showed at least some "curvature" in the plot indicating compressible peat layers, but the other 3 holes produced "straight lines" indicating sites where the soil material was apparently more rigid.The tests producing "straight lines" (Table 1: BH7, BH8 and BH10) were only interpreted using the steady-state methods.
The tests producing "curved lines" have been interpreted using the transient KGS and Dax methods, but for comparative reasons have also been analyzed using the steady-state van Beers and Hvorslev methods.The computed K h values fall within the range of 0.10 m d −1 to 1.69 m d −1 which are normal values for peat land areas in The Netherlands (Weerts, 1996).The consistently lower values obtained with the KGS approach in relation to the Dax method may be due to the difficulty in determining accurate S values for the latter method.
A comparison can be made between the obtained K h values, the compressibility at the site and the decomposition degree of the peat.The latter can be expressed using the von Post scale.The boreholes with higher conductivities and a higher compressibility interpreted with the transient KGS and Dax methods have a higher decomposition degree (Table 1: BH 5,6,9,11,12).Typically, these slug tests were carried out in well decomposed peat which has a muddy and plastic appearance that also contained mineral components.The holes with lower conductivities and a more rigid soil analyzed with the steady-state van Beers and Hvorslev formulae have a relatively low decomposition degree (BH 7, 8, 10).These tests were carried out in more coherent soils made up of clay, silt and peat or where fresh, largely intact peat was found.
The inappropriate use of steady-state formula to interpret slug tests in compressible peat has wrongly led to the belief that Darcy's law is not valid for these types of materials, but on the other hand full proof has also not been given that this most important law for groundwater can be applied in peaty environments.In other words, is the hydraulic conductivity K really independent of the groundwater head differences in the groundwater flow system?Waine et al. (1985) presented evidence of a higher hydraulic conductivity K of well-decomposed peats with an increase in head gradient for steady-state experiments set up in columns in the laboratory.The way in which Waine did the experiments also meant that an increase in gradient increased the mean head in the column.Nevertheless, they found that the effect of the increase in gradient was rather small.For a three-fold increase in head gradient, the hydraulic conductivity only increased by about 30 %.The slug tests carried out for the current study were conducted for different initial groundwater heads.No trend was detected between the hydraulic conductivity and the magnitude of the heads applied.Based on this information the statement can be made that the K may not be entirely independent of the groundwater head differences in a system.However, in a natural peaty environment the head differences are generally small and will hardly affect the K, justifying the application of Darcy's law in analytical and numerical model computations.In modeling on the other hand, even though the K is affected, there are usually larger errors arising from uncertainties in the hydro-geological conceptualization and model input data.

Modeling results
The groundwater modeling at the research area incorporates the heterogeneities of the Complex Confining Layer (CCL) that were disclosed from the fieldwork.The modeling results have been compared with the outcome of a simplified groundwater model whereby the CCL is represented as a single homogeneous layer.In this way the effect of heterogeneity on modeling results has been established.For the comparison, groundwater models have been constructed which have been calibrated using the available groundwater heads (levels) at observation wells.However, using only heads for model calibration may lead to solutions that are nonunique (Konikow and Bredehoeft, 1992).To obtain a correct solution, thorough attention has been paid to the hydrogeological conceptualization and data input for the models.

Layer schematization for conceptual model
The heterogeneities in the CCL made it clear that a one layer schematization of the CCL is too simplified, and at least a distinction has to be made between the muddy and peaty layers and the more clayey layers.Based on the cross section obtained from the drilling activities (Fig. 5), the schematization of the CCL in the model area included a representation by 4 model layers (Fig. 8).The alteration of sandy aquifers and clayey aquitards below the CCL was represented by another 9 model layers, bringing the total number of layers to 13 (Fig. 9).The base of the aquifer system, which is equivalent to the base of model layer 13, is located at a depth of about 310 m below surface in the so called Maassluis Formation.
The option to limit the model to the CCL was considered but discarded, since groundwater head data that could be used to set up the base of the model in the first sandy aquifer below the CCL were thought to be too inaccurate, in particular for the steady state model.Such levels would be a result of a Kriging exercise based on groundwater head measurements taken at boreholes outside the model area (Fig. 1).In addition, the model layer schematization and data input below the CCL could be based on the information contained in an existing well-calibrated regional model (Cheng, 2004).

Boundary conditions for conceptual model
The groundwater model area was selected within the research area and also incorporated the experimental field that measures 50 m by 300 m (Fig. 1).The groundwater headcontrolled conditions had to be assumed for all the models layers in the lateral boundaries of the model area (which had an area of 700 m by 900 m).Modeling an area larger than the experimental field where most of the data were collected was allowed as can be deduced from lithological considerations.The lithological information at the experimental field revealed a heterogeneous lithology in terms of layer thicknesses and even layer discontinuity, but -apart from the  top layer -there was also a consistent alteration of a more peaty upper part, a clayey middle section, and a peaty bottom section of the CCL.This latter structure was also visible in most of the boreholes in the vicinity of the experimental field (Fig. 1).Therefore, the model of 700 by 900 m accurately represents the lithological character of this peat area in The Netherlands as revealed by the detailed descriptions at the experimental field.

Groundwater balance
A groundwater balance can be defined for the entire groundwater model area including all model layers.The focus of attention in this paper is on the CCL; therefore, the balance for this confining layer has been introduced as follows: where Q prec (m 3 d −1 ) denotes the recharge at the groundwater table in the CCL as a result of precipitation, Q cap (m 3 d −1 ) expresses the capillary rise from the table, Q surfin and Q surfout (m 3 d −1 ) refer to the groundwater inflow and the outflow at the ditches, respectively, and the Q up and Q down (m 3 d −1 ) are the vertical upward and downward flow exchanges between the CCL and the underlying Pleistocene sandy aquifer, respectively.The term V CCL / t (m 3 d −1 ) describes the change in water storage per time step t in the CCL.The groundwater model and its balance refer to the part of the CCL that is saturated, which is by far the largest part of this unit.In addition, a model has been built for the unsaturated top part of the CCL, allowing the assessment of the precipitation surplus (Q prec -Q cap ) required for the saturated model (see Sect. 5.4.3).With regard to the direction of groundwater movement at the research area, the flow through the CCL is downward, which means that the Q up can be neglected.A thorough analysis of phreatic groundwater heads in the CCL and confined heads for the sandy aquifer indicated this flow direction (Fig. 3).

Model code and data input
The selected software for the modeling at the research area was the well-known code PMWIN-MODFLOW 5.3 which is based on Darcy's law.The Recharge and River packages as part of this code were implemented.The grid for the MOD-FLOW model was designed with a constant cell size of 5 m by 5 m.This relatively small cell size relates to the small widths of about 50 m of the plots of grass land characterizing the research area (Fig. 2).In order to compute and calibrate the groundwater heads in the plots bounded by ditches with sufficient accuracy, this small cell resolution was required.
Since very detailed field data were not available, the selection of an even smaller grid size than 5 m by 5 m was not necessary.

Hydro-geological parameters
A Digital Elevation Model (DEM) with a resolution of 5 m by 5 m provided the information to input into the model the elevations of the land surface that varied between −1.9 m and −2.7 m.Based on fieldwork data, the elevations of the top and bottom of the model layers making up the CCL were inserted into the model.The bottom of the CCL varies between −8.4 m and −9.1 m.Borehole data of the database residing with the Institute of Applied Sciences (TNO) were interpolated to obtain adequate elevations of the model layers, corresponding with the sandy aquifers and clayey aquitards below the CCL (Cheng, 2004).
Slug test experiments carried out during the fieldwork, permeameter laboratory tests completed by the Wageningen University and Research Centre, various literature that discuss the permeability of the CCL (e.g.Weerts, 1996), and pumping tests completed in nearby pumping stations for domestic water supply were considered in providing the data to set up the hydraulic conductivities and storage parameters for the different model layers (Table 2).For example, the horizontal hydraulic conductivities ranged from 0.0024 m d −1 for clay layers to 1.3 m d −1 for peaty units.To satisfy anisotropic conditions, ratios between the horizontal and vertical hydraulic conductivities, varying from 1 to 2.7, were adopted for the peaty and clayey model layers in the CCL (Weerts, 1996).

Time and model boundary heads
Steady-state and transient models were prepared for specific periods.The steady-state model runs on the input of data for the year 1990.This year reflects average meteorological conditions and the model simulates average groundwater heads and flows (Cheng, 2004).The transient model was built for the winter period 2009 to 2010 adopting a stress period of 10 days.Groundwater head data collected at the experimental field were available to calibrate this model.The model was also extended to simulate seasonally varying flow in the CCL from 2006 to 2010 and from 2002 to 2010 using excellent sets of historical data.The extension was also necessary in order to visualize the complete pathlines of water particles through the CCL and to be able to compute corresponding travel times.
Field-based measurements defined the boundary and initial groundwater heads for the modeling area.For the peaty layers with a phreatic response in the CCL the heads, ranging from −2.3 m to −2.9 m, were obtained from so called GxG maps.The phreatic surface, which is rather irregular in the ditch and land plot landscape of the model area, was taken from these groundwater head maps which have a grid resolution of 25 m by 25 m.The GxG maps are official maps compiled on the basis of groundwater head (level) measurements in boreholes monitoring the phreatic layer, followed by the application of interpolation techniques.For the Pleistocene sandy aquifer and the deeper aquifers, the TNO database provided the groundwater heads for the model, which ranged from −2.9 m to −3.9 m.This information needed to set up the boundary and initial groundwater heads for the deeper model layers was estimated by interpolating the data available for the few boreholes present in the area.

Groundwater recharge and the surface water system
An unsaturated model for the root zone, implemented in an Excel spreadsheet, supplied the input data for groundwater recharge from precipitation.This model, considering as main input parameters soil characteristics, precipitation, and potential evapotranspiration, computes net recharge as the balance between downward recharge and capillary rise.Precipitation data for the rainfall station within the research area itself (Fig. 1) and evapotranspiration data of the De Bilt weather station were input into the unsaturated zone model and subsequently the corresponding net recharges were computed.For example, for the observation period selected for the steady model, the computed net recharge varied from a high of 2.3 mm d −1 to a low of −1.7 mm d −1 , yielding a realistic average recharge value of 0.41 mm d −1 as the input for the model.The Water Board HDSR provided surface water information including data relating to the River Oude Meije and the extensive ditch system.The data comprise river and ditch characteristics like bed hydraulic conductances, open water levels, and elevations of the bottom of these surface water courses.The open water levels at the river and ditches in the model area are controlled by HDSR and maintained at a fixed position for winter and summer conditions.

Sensitivity analysis and model calibration
The data input for the models was followed by a sensitivity analysis aiming to identify the parameters which have a large effect on modeling results including the computed phreatic Hydrol.Earth Syst.Sci., 15, 3017-3031, 2011 www.hydrol-earth-syst-sci.net/15/3017/2011/  and confined groundwater heads.The analysis showed that increases or decreases of 50 % in the values for the hydraulic conductivity of the CCL and the groundwater recharge resulted in phreatic groundwater head changes of less than 0.1 m.Confined groundwater head changes were even below 0.01 m.The hydraulic conductivity is a bit less sensitive than the recharge.The transient models require the input of storage parameters like the specific yield.The value of this parameter, which increased 50 %, consequently resulted in a phreatic head change of 0.08 m.This indicates a sensitivity of a similar order to the sensitivity of the conductivity or the recharge.Referring to the numbers mentioned the sensitivity of model parameters is not that large.This could be attributed to the control which the open water levels at the ditches exert on the groundwater heads.The models were calibrated taking into account the results of the sensitivity analysis, which meant that parameter values were changed in order to try to minimize the difference between the groundwater heads computed with the model and the groundwater heads measured in the field.Although the hydraulic conductivities are less sensitive than groundwater recharge, the main model adjustments concerned the former parameter.The reason is that the values for the conductivities are believed to be less certain than the values for the recharge.
During the actual calibration work with the steady state model, the phreatic groundwater heads computed on the bases of a particular set of parameter values were compared with the field-based phreatic heads shown on the GxG maps.For the transient model, computed phreatic and confined groundwater heads could be compared with groundwater heads measured during the fieldwork in the upper peaty layers of the CCL and in the sandy aquifer.Model calibration resulted for final mean absolute differences in computed and measured groundwater heads of less than 0.05 m, in particular in upgraded model conductivities ranging from 0.0001 m d −1 for the clayey parts of the CCL up to 1.3 m d −1 for the peaty layers (Table 3).

Comparison of detailed and simplified models
Results of the calibrated detailed models with a heterogeneous CCL were compared with the outcomes of simplified models whereby the CCL is represented as a single homogeneous layer.To obtain a representative simplified model the upper four layers of the detailed model were merged into a single layer with equivalent values for parameters including the hydraulic conductivities and the effective porosity.The equivalent horizontal conductivity was computed as the total horizontal transmissivity of the CCL divided by the thickness of this layer, resulting in a K h value of 0.961 m d −1 .The equivalent vertical conductivity was elaborated from the thickness of the CCL divided by the total vertical resistance across this unit.Some aerial variation in resistance had to be taken into account, leading to the computation of equivalent values with magnitudes of K v1 = 0.0114 m d −1 , K v2 = 0.0061 m d −1 and K v3 = 0.0013 m d −1 .For the equivalent porosity an average of the 4 model layers of 0.26 was considered.The assignment of equivalent values to the simplified model meant that the unique effect of heterogeneity could properly be determined when comparing model results.
Model results that can be compared comprise pathlines and travel times.For either heterogeneous or homogeneous conditions, pathline plots demonstrate the different trajectories that water particles may have followed through the CCL; and travel time computations indicate the time it takes a water particle to flow through the CCL to the sandy aquifer or the ditch system.The prepared detailed and simplified models for the research area are representative for peat areas where groundwater flow is primarily directed in a downward direction.Additional hypothetical models have been prepared to study the pathlines and travel times for typical upward flow through the CCL.These latter models have been prepared given that in other peaty environments in the western part of The Netherlands, the dominant flow of the water particles is upward (Gonzales et al., 2009).

Steady state model
Steady state models for the research area were used to obtain pathlines and travel times for the typical case of downward flow through the CCL.The models assume that average hydrological conditions are maintained over a long time.The model results show that for both the heterogeneous and homogeneous cases, the flowpath of water particles are similar and the direction is almost vertically downward (Figs. 10 and  11).When the particles arrive in the sandy aquifer, the flow becomes more horizontal.For the heterogeneous case the travel time of water through the saturated CCL was approximately 7 yr, and for the homogeneous case a water particle was underway in the CCL for 9 yr.

Transient model
Transient models that have been engaged to simulate pathlines and compute travel times did not assume the average state, but they did take into account the seasonally varying hydrological conditions.Therefore, they allowed a more  the travel time through the saturated CCL is on the order of 4 yr, but for the latter case the water particles have not even reached half way through the covering layer in this 4yr period.The travel times for the transient models are also less than for the steady models, indicating that the differences in head distributions between the models play a prominent role.These differences can be attributed to the average hydrological conditions and constant water levels assumed for the steady simulation as well as the seasonally varying conditions and levels adopted for the transient state.
How can the longer travel time for the homogeneous transient model be explained?In the homogeneous model the hydraulic conductivities of the different layers are replaced by a single representative conductivity.This means that the top layer of clayey and silty peat obtains a much higher conductivity, resulting in an increase in groundwater flow through this layer towards the ditches.The flow from the top (model) layer to the second layer will subsequently be much less than for the heterogeneous case, causing smaller flow velocities at the interface between these layers.The MODFLOW associated code MODPATH which was engaged to determine pathlines and travel times, uses the velocities at the interfaces to compute the velocity distribution within the layers.This meant that in the whole second model layer, the velocities of the homogeneous model are small (Pollock, 1988).The small velocities in the second, and also the third, model layer lead to the computation of long travel times for the homogeneous case.

Typical upward flow through the CCL
Transient hypothetical models have been made to generate pathline patterns and compute travel times for upward flow through the CCL.The upward flow models have the same set-up, boundary, and initial conditions as the downward flow models, the only difference being that the groundwater heads in the confined aquifer are higher than the heads in the phreatic aquifer.For comparative reasons, the absolute differences in confined and phreatic groundwater heads for these models were similar than for the models simulating downward flow.To visualize the flow through the CCL, water particles had to be placed both at the water table and in the Pleistocene sandy aquifer.The model results indicate that a complex flow pattern develops, culminating in a prominent outflow at the ditches (Figs. 14 and 15).Although upward flow is dominant, also downward flow exists in the upper parts of the CCL.The heterogeneous case also clearly shows in the second peaty model layer and centrally between the ditches, the development of an area with nearly stagnant groundwater flow.The travel times calculated for upward flow are comparable with those for downward flow.For the heterogeneous case the travel time for a water particle to move from the sandy aquifer into the ditch is on the order of 2 to 4 yr, and when homogeneous conditions are considered the travel time to the ditch is longer -ranging from 5 to 8 yr.Not surprisingly, the travel times are shortest for particles that move upward in the CCL below the ditches.Particles that flow upward centrally between the ditches where the stagnant area is located take the longest time to finally reach the surface water system.

Conclusions
The shallow surface water and groundwater system in a typical peaty lowland area in The Netherlands could successfully be analyzed using a combination of fieldwork and modeling work.Uncertainties in modeling arising from the heterogeneity in peat soils and the applicability of Darcy's law could be eliminated through a proper model set up.The model itself provided useful information on the flow path patterns and travel times of groundwater through the Complex Confining Layer (CCL) as well as the flow exchange with the ditch system.Our fieldwork and previous investigations (e.g.Weerts, 1996) indicated that the CCL indeed deserves its reputation of being a complex layer.Borehole drilling revealed that in the vertical and horizontal direction the CCL is heterogeneous.In the vertical direction, layers of peaty and silty clay, mud and peat, clay, and again peat succeed each other in a downward direction at the research area.The CCL rests on top of fine eolian sand which grades into coarser fluvial sands at lower levels.In the horizontal direction the various layers are continuous in many places, but also tend to fade out, offering additional complexity for groundwater flow.
Field slug tests indicated that an apparent variable K h was observed for humified peat soils, but this is thought to be due to the inappropriate use of steady-state formula for transient experiments.Even though these expressions are normally used for transient slug tests, they were found to be invalid for these kind of peat soils.During the tests the humified peat soils introduced a large transient effect as a result of the compressibility of the peat and its plastic nature.Based on interpreted elastic storativity values of the peat, varying from 0.0001 to 0.001, the K h of the humified peat soil could be computed with an equation for transient conditions and a semi-analytical model yielding values in the range of 0.24 m d −1 and 1.69 m d −1 .On the other hand, the slug tests gave no clear evidence of the dependency of the hydraulic conductivity with the applied hydraulic heads.In combination with earlier results from laboratory tests (Waine et al., 1985), the evidence proved that normal (pressure) head differences in a natural groundwater system hardly affect the hydraulic conductivity.Therefore, Darcy's law can be used for the typical peat soils as tested in The Netherlands.
The modeling work focused on the comparison of a detailed model for the peaty CCL and a model with lumped parameters.The detailed model takes into account the heterogeneity of the CCL, which is represented by 4 model layers; but the lumped model considers the CCL as an apparent homogeneous entity.The transient models running for simulation periods from 2006 to 2010 and 2002 to 2010 proved to be adequate to visualize the pathline patterns in the CCL for downward flow in the research area and for upward flow encountered in other similar areas in the western part of The Netherlands.As a result of the similar hydrological conditions adopted at the boundaries of the CCL, the flowline patterns for the model representing the heterogeneous case are quite similar to those generated by the model with a homogeneous CCL.Flow patterns made with the heterogeneous model tend to be more accurate when looking at flow details, for example when inspecting the flow exchange between the CCL and the surface water system.
An interesting result from the modeling exercises are the indications of the travel times of water particles when flowing through the peaty CCL.For the heterogeneity case where the CCL is taken into account, the model computes travel times of water particles through the CCL which are considerably shorter than for the case whereby the Complex Confining Layer is represented by a single homogeneous model layer.Travel times downward or upward through the CCL for the heterogeneous model range between 2 to 4 yr, whereas the model with homogeneity tends to compute times on the order of 5 to 8 yr.The large over estimation of the travel time for the homogeneous case applies to the model area as part of the research area, but elsewhere in the HDSR control area the discrepancies in travel time between homogenous and heterogeneous cases will be be different.The conclusion is that groundwater models that are based on the representation of the CCL with one homogeneous model layer are not suitable for assessments on groundwater transport where travel times play an important role.Taking one homogeneous layer into account could result in the computation of travel times that are even double the values for the heterogeneous case.In particular when models are considered for the simulation of contaminant transport, models with a homogeneous CCL should not be used.

Fig. 1 .
Fig. 1.Location of the research area in The Netherlands (left) within the HDRS control area (top right).(a) Research area: w = Wageningen experimental area; e = current experimental area; r = rainfall station; m = model area; polygons (e.g.1265) are TNO boreholes with lithological description; squares (e.g.0126) are TNO boreholes with head measurements.

Fig. 2 .
Fig. 2. Location of the installed equipment and experiments carried out during fieldwork.
Fig.3Measured groundwater levels in three piezometers (pressure transducer PT 2 measures confined levels and PT 3 and PT 4 measure phreatic levels.PT 1 had to be discarded).

Fig. 4
Fig. 4 Subsoil description of the deeper boreholes

Fig. 3 .
Fig. 3. Measured groundwater levels in three piezometers (pressure transducer PT 2 measures confined levels and PT 3 and PT 4 measure phreatic levels.PT 1 had to be discarded).

Fig. 3
Fig. 3 Measured groundwater levels in three piezometers essure transducer PT 2 measures confined levels and PT 3 and PT 4 measure phreatic levels.PT 1 had to discarded).

Fig. 4
Fig. 4 Subsoil description of the deeper boreholes

d
Fig. 6 Geophysical investigation and interpretation

Fig. 5 .
Fig. 5. Schematization of the CCL based on borehole descriptions.The average groundwater level for an average year is at about 40 cm below the ground level.

Fig. 6 .
Fig. 6.(a) Field resistivity plot with 3-layer interpretation (given by the shape of the obtained curve).(b) Geo-electrical profile.Geophysical investigation and interpretation.

Fig. 10 Fig. 10 .
Fig. 10 Steady pathlines for downward flow for the heterogeneous case (b) lines are pathlines; each dot represents one year (b) blue boxes are ditches (c) blue arrows are velocity vectors indicating direction and magnitude of water movement (d) note that velocity does not decrease in the clay (layer 3) due to large groundwater level difference across this layer

Fig. 11 .
Fig. 11 Steady pathlines for downward flow for the homogeneous case (a) lines are pathlines; each dot represents one year (b) blue boxes are ditches (c) blue arrows are ve vectors indicating direction and magnitude of water movement

Fig. 12 Fig. 12 .Fig. 13 .
Fig. 12 Transient pathlines for downward flow for the heterogeneous case (a) flow simulation 2006 to 2010, each dot represents one year (b) blue boxes are ditches (c) blue velocity vectors indicating direction and magnitude of water movement

3Fig. 14 Fig. 14 .
Fig. 14 Transient pathlines for upward flow for the heterogeneous case (a) flow simulation 2002 to 2010, each dot represents 4 months (b) blue boxes are ditches (c) blue arrow velocity vectors indicating direction and magnitude of water movement

)Fig. 15 .
Fig 15.Transient pathlines for upward flow for the homogeneous case ) flow simulation 2002 to 2010, each dot represents 4 months (b) blue boxes are ditches (c) blue arrows are velocity vectors indicating direction and magnitude of water movement

Table 1 .
Horizontal hydraulic conductivity K h (m d −1 ) computed from the slug tests.

Table 2 .
Horizontal and vertical hydraulic conductivity from different sources.

Table 3 .
Horizontal and vertical hydraulic conductivities for the model layers representing the CCL.