Response to recharge variation of thin lenses and their mixing zone with underlying saline groundwater

S. Eeman, S. E. A. T. M. van der Zee, A. Leijnse, P. G. B. de Louw, and C. Maas Wageningen University, Environmental Sciences Group, Soil Physics, Ecohydrology and Groundwater Management, P.O. Box 47, 6700 AA, Wageningen, The Netherlands Deltares, Dept. Soil and Groundwater, P.O. Box 85467, 3508 AL, Utrecht, The Netherlands KWR Watercycle Research Institute, P.O. Box 1072, 3430 BB, Nieuwegein, The Netherlands

by these variations in recharge, contrary to the mixing zone thickness.The mixing zone thickness is compared to that of a Fickian mixing regime.A simple relation between the travelled distance of the center of the mixing zone position due to variations in recharge and the mixing zone thickness is shown to be valid for both a sinusoidal recharge variation and actual records of daily recharge data.Starting from a step response function, convolution can be used to determine the effect of variable recharge in time.For a sinusoidal curve, we can determine delay of lens movement compared to the recharge curve as well as the lens amplitude, derived from the convolution integral.Together the proposed equations provide us with a first order approximation of lens characteristics using basic lens and recharge parameters without the use of numerical models.This enables the assessment of the vulnerability of any thin fresh water lens on saline, upward seeping groundwater to salinity stress in the root zone.

Introduction
Rain-fed areas may suffer from salinity in the root zone when salt groundwater is found at shallow depths.In low-lying coastal zones, saline water is often present at shallow depth, due to a history of flooding (Vos and Zeiler, 2008;Post, 2004), marine transgressions and sea spray (Stuyfzand and Stuurman, 1994).In such areas, infiltrating rain water is the only source of fresh water, forming and maintaining thin fresh water lenses floating on top of the saline groundwater.Infiltration of rainwater is limited by upward seepage of the saline groundwater when the soil surface is below sea level as found in deltaic areas like the Netherlands (De Louw et al., 2011;Maas, 2007;Oude Essink et al., 2010), and foreseen due to future relative sea level rise in for example the deltas of the Nile and Mississippi (Jelgersma, 1996).The surface area of such coastal zones and delta regions is considerable and they are generally densely populated.Therefore, and in view of the usually good water supply and soil fertility conditions, their agricultural and ecological importance is significant.The integrity of fresh water lenses is threatened by both the expected sea level rise (Day et al., 1995;Lebbe et al., 2008) and the drainage of soil for agricultural reasons, which have caused subsidence of these low-lying areas and therefore increased the saline ground water pressure ( Van de Ven, 2003).Studies by for example Zeidler (1997), Day et al. (1995), Habibullah et al. (1999), andOsterkamp et al. (2001) conclude that salinization in lowlying coastal areas in different parts of the world will cause changes in existing ecosystems and in some places loss of arable land over the next decades if no measures are taken.In all cases, water management will play a dominant role.To minimize adverse effects of salinization, water management and land use planning should be based on a proper understanding of the behavior of thin fresh water bodies on saline groundwater and their response to changes in recharge, whether seasonal, longer term, natural or human induced.
Related to the behavior of thin lenses is the problem of the classical dune system, a thick lens of fresh water floating on top of more or less stagnant saline ground water.Fresh water bodies in dunes have received a lot of attention from Badon Ghyben (1888) and Herzberg (1901) and, more recently, by many others who derived analytical solutions for sharp interface steady state situations (e.g.Bakker, 2000;Bruggeman, 1999;Fetter, 1972, Van der Veer, 1977).Hantush (1968), Van der Veer (1977) and Bruggeman (1999) also provided sharp interface solutions for moving fresh water lenses in saline aquifers under different conditions.Maas (2007) derived a steady state approximation for lenses on top of upward flowing saline water, which is very similar to the relatively thin lenses considered in this paper.
For relatively thin lenses, a sharp interface approach is not appropriate since the mixing zone between the fresh and saline water is thick compared to the total lens thickness, as was demonstrated by De Louw et al. (2011).This mixing zone between fresh and saline water has been studied for several steady state situations by for example Henry (1964), De Josselin de Jong and Van Duijn (1986), Paster and Dagan (2007), and Abarca et al. (2007).
Many case studies use numerical simulations to investigate mixing zones at different spatial and temporal scales.Of main interest is to identify the dominating process(es) that influence the mixing zone and the delay of the response of the mixing zone, but results appear very dependent on the environmental conditions.Underwood et al. (1992) found for atolls that mixing was controlled by short-term fluctuations, driven by tidal pulses, whereas recharge drives the long term average ground water flow that determines lens thickness.However, Kiro et al. (2008) concluded that the lowering of the Dead Sea caused a delayed reaction in the mixing zone.This movement, very slow and small compared to tidal fluctuations, apparently needs more time to affect the mixing zone.In agreement with Kiro et al. (2008), Andersen et al. (1988) saw that small but long-term changes in groundwater level in a drinking water well area can significantly influence the position of the mixing zone, whereas annual cyclic changes did not have any influence.Cartwright et al. (2004) investigated the response of the mixing zone to tides and waves on a beach.Contrary to Underwood et al. (1992), they found that the effect of the diurnal tidal pulse (0.5 to 2 m sea level rise) did not influence the Introduction

Conclusions References
Tables Figures

Back Close
Full position of the center of the mixing zone.However, they observed a nearly immediate and strong effect of the wave-induced groundwater pulse on the average position of the mixing zone during moderate storm events (wave heights ∼4.5 m during 1-3 days).The higher amplitude and lower frequency enhance each other's effect on increasing the mixing zone.Cartwright et al. investigated much smaller spatial (horizontal and vertical) and temporal scales than Underwood et al. (1992) and Kiro et al. (2008), and neglected thickness changes of the mixing zone.As these illustrations indicate, the response depends on the physical system parameters, the spatial scale and the duration and amplitude of temporal variations of the flux or water level.It is plausible, that it also matters whether such variations are forced at the perimeter of the system (e.g.tidal variations), and may therefore attenuate with distance from this perimeter, or are forced along the entire top boundary.The latter is of interest, if one wishes to consider variations in time of precipitation that lead to fresh water recharge.
Whereas fresh water lenses in larger dune areas generally have reaction times of years to even tens of years (e.g.Vaeret et al., 2011;Oude Essink, 1996), for thin lens systems, we expect fast responses of the position of the mixing zone to recharge changes.The main responses of interest are the thickness of the lens and of its mixing zone.For a lens changing in thickness, flow lines are predominantly vertically oriented, perpendicular to the mean interface (Eeman et al., 2011).This agrees with observations (De Louw et al., 2011) and implies that mixing is controlled by longitudinal dispersion, which creates a thicker mixing zone compared to a transversal dispersion/diffusion dominated stable mixing zone.The combination of the position and width of the mixing zone determines together with the capillary characteristics of th soil, whether saline water reaches the root zone during the growing season, and concern about the impact of root zone salinity on primary production is one of the motivations of this research.
Our objective is to relate changes of fresh water lenses and their mixing zone to hydrological characteristics of a field site and recharge variation.This relation can be Introduction

Conclusions References
Tables Figures

Back Close
Full used to determine under which circumstances water quality is affected in the root zone.
Lens properties that dominantly influence the depth at which saline water is found, are the lens thickness and the vertical extent of its mixing zone.We analyze the effect of recharge variations with different duration and intensity on lenses with a different thickness and mixing zone thickness, using numerical simulations.

Flow domain and system equations
The modeled lens system is comparable to that in our previous work (Eeman et al., 2011).It is a vertical cross section of the flow domain between two ditches or drains Introduction

Conclusions References
Tables Figures

Back Close
Full (Fig. 1a).Recharge with fresh water (uncolored) occurs at the upper boundary, while vertical upward seepage of saline water (grey-colored) occurs at the lower boundary.All water entering the flow domain leaves through the ditches or drains.Such situations are common in deltaic lowlands, and field observations (Fig. 1b) illustrate that the resistivity, a measure of the salinity of the pore water, decreases gradually with increasing depth.The resistivity sequence found can be completely attributed to the salt concentration of the saline groundwater, which is for this case around two thirds of the salinity of sea water.
In our analysis, the volume of fresh water in the lens includes the fresh water in the unsaturated zone.This is a reasonable approach because the unsaturated zone is thin (generally around 1 m), but also in view of our interest in the fresh water availability for plants.In the text, we refrain from mentioning this issue all the time, for brevity.
We briefly give the equations for the water-saturated domain, and for more detail, we refer to our earlier work (Eeman et al., 2011).Symbols and units are in the nomenclature.The mass balances for the fluid and salt are, respectively, given by where ρω is the mass density of the salt in the fluid phase and ρωq i are the convective components of the salt flux in the i -direction (i = x,z).J is the dispersive flux as defined by e.g.Bear (1972).The system is assumed to be incompressible.The equation of state gives the fluid density ρ as a linear function of the salt mass fraction ω (Weast, 1982).
We assume that the porous medium is isotropic and take the z-coordinate positive in the downward direction.The components q x and q z of the specific discharge q of the Introduction

Conclusions References
Tables Figures

Back Close
Full fluid are given by Darcy's law where κ and µ are constants, i.e. we consider a homogeneous medium.A constant viscosity µ is justified as it changes by less than 2 % between fresh water and seawater (Weast, 1982).
Based on all system equations, we identified 8 dimensionless groups (Table 1) that define the steady state system in Eeman et al. (2011).Lens thickness is particularly sensitive to changes in mass flux ratio M and Rayleigh number R (as is shown in the figure attached to Table 1).The mass flux ratio M is the ratio between the mass fluxes of the lower and upper boundary, the Rayleigh number R is the ratio of the characteristic density induced flux and the reference specific discharge P (Table 1).
To assess the response of lens and transition zone to variations in recharge, we adapted the top recharge boundary condition used by Eeman et al. (2011).We use sinusoidal recharge variations (P s ) to model variations in recharge, while the average P replaces the constant P in M and R.This enables us to simply vary the main features of the recharge, such as average, duration, and intensity.Moreover, the absence of abrupt changes in sinusoids is numerically more attractive.
The boundary condition for the total mass flux at the upper boundary (fresh water precipitation) is given by We obtain the following two additional dimensionless groups by combining Eq. ( 5) and the reference values as elaborated in the steady state analysis of Eeman et al. ( 2011) and

Conclusions References
Tables Figures

Back Close
Full Where P is the average of P s , and the groups represent the dimensionless frequency and amplitude, respectively.In our analysis we will focus on these two groups, and the mass flux ratio M and the Rayleigh number R.

Model schematization and parameter values
Only one half of a lens system is considered for reasons of symmetry.Besides the precipitation at the upper boundary, the boundary conditions that we used are shown in Fig. 2: a (saline) seepage flux at the bottom boundary, a pressure boundary along the ditch and closed vertical left and right boundaries, discussed in detail by Eeman et al. (2011).SUTRA (Voss and Provost, 2008) was used to carry out density-dependent numerical simulations.Spatial discretization of the quadrilateral elements was chosen such that the numerical dispersion was much smaller than the mechanical dispersion for all simulations.The element size in the vertical direction was 0.1 m for the upper 1.5 m of the domain (Fig. 2, zones 1 and 3) and 0.2 m for depths larger than 1.5 m (Fig. 2, zones 2 and 4).The element size in the horizontal direction was 0.1 m in the outflow area near the ditch (Fig. 2, zones 1 and 2).In the infiltration area horizontal element sizes of 0.4 m were used (Fig. 2, zones 3 and 4).
Temporal discretization was controlled by criteria that limit the allowable changes in pressure, saturation and concentration per time step, and time step sizes were adapted accordingly.An overview of reference parameter values is given in Table 2.For the initial conditions, steady state lenses were used that were obtained with a constant recharge rate of 1 mm day −1 .
The volume and thickness of rainwater lenses and the thickness of the mixing zone were quantified using spatial moments (Eeman et al., 2011).Spatial moments efficiently summarize the numerical results and have been used widely in other contexts, e.g.Acharya et al. (2005) and Paster and Dagan (2007).Govindaraju and Das (2007) discuss formulation and use of spatial moments.From the zeroth spatial moment of the salt mass fraction over the modeled domain, the total volume of fresh water in the lens and transition zone can be inferred.The spatial moments of the vertical derivative 1445 Introduction

Conclusions References
Tables Figures

Back Close
Full of the salt mass fraction (inset Fig. 1a) provide information on the position, thickness, and shape of the mixing zone and can also be used for field measurements, as shown by De Louw et al. (2011).
To assess the response of the lenses under different conditions to net recharge, we varied the lens characteristics and the recharge.We varied average lens thickness by varying M and R (Tables 1 and 3), which leaves all other groups constant (using reference parameter values as given in Table 2).Only C, the dimensionless conductance of the soil-ditch interface, changes as we change R. Since the influence of C is very small (inset Table 1), its variation has a negligible impact on the results.We focus on thin lenses because these are relevant for primary production in relation to salinity stress.For lenses thicker than 3 m and commonly observed recharge regimes, water with significant amounts of salt is very unlikely to reach the root zone (Eeman et al., 2011).For lenses that are thinner than about 0.8 m, the salt stress is likely to be so severe, that it limits plant growth to halophytic species.
To parameterize the sinusoidal recharge P s we used variations from a reference P s fitted to average Dutch weather data (De Bilt, 1971-2000, provided by the Royal Netherlands Meteorological Institute KNMI), see Fig. 3a and the bold values in Table 3.
To quantify the influence of the amplitude and frequency of P s on lens response, the amplitude was varied between 0.23 and 23.0 mm day −1 , while the period was varied from 1 yr to 1 week (Table 3).P , defined as the average of P s was kept at 1 mm day −1 .
The strongly fluctuating weather data from the KNMI (Fig. 3) are used to validate the approximation of the thickness of the mixing zone that we will propose (Sect.3.2).

Mixing
Particularly if the thickness of the mixing zone between fresh and salt water is of concern, the choice of the dispersivities is important.Previously (Eeman et al., 2011), we used dispersivities of 0.25 (longitudinal) and 0.05 (transversal), which comply with reported macro dispersivities for aquifers (Gelhar, 1992;Kaleris, 2002Kaleris, , 2006)).The value of the dispersivities that follows from Fiori and Dagan (2000) would easily be Introduction

Conclusions References
Tables Figures

Back Close
Full one (or more) order(s) of magnitude smaller.The macroscopic values represent pore scale mixing, but also larger scale spatial variability of soil that does not necessarily lead to true mixing.Mostly, those values are derived for two scales of variability (pore scale and scale of variability of hydraulic conductivity in porous formations).Due to exchange of solute between stream tubes with different velocities, stream tube interfaces enhance true mixing (Janssen et al., 2006), similar as mobile/immobile exchange (Van Genuchten and Dalton, 1986;Parker and Vallocchi, 1986;Bolt, 1982).Larger dispersivities, even after considering their reliability (Gelhar et al., 1992), may be caused by model (including dimensionality), fitting, and experimental bias.For instance, the mixing may well occur in the extraction wells: if wells have a relatively large screen, water from different strata is mixed.Non-invasive techniques, as used for the data of Fig. 1b, may lead to apparent mixing by averaging through a limited spatial support of the technique.For instance, in Fig. 3b, we compare simulated and measured mixing zones for two sites at Schouwen-Duiveland (South West Netherlands) that we monitored for the past two years (De Louw et al., 2011).Figure 4 illustrates that a dispersivity α l of 0.25 m (for a moving interface) reflects the observed mixing zone thickness of 1 to 2 m better than the smaller α l of 0.05 m.However, these larger dispersivities are obtained by fitting a homogeneous 2-D model that disregards the actual layering at the sites.
A relatively large effective dispersivity could be caused by different processes.For non-uniform media, the presence of mobile and immobile water regions affects mixing.This leads to a larger effective dispersion coefficient of the form (Parker and Valocchi, 1986)

System delay using convolution
To investigate the delay between the changes in recharge and thickness response of thin lenses, we adopted a standard systems analysis, using the impulse response function and convolution of this function with the recharge signal (e.g.Olsthoorn, 2008).Practical solutions using convolution can for example be found in Maas (1994) and Bruggeman (1999).Although this approach is only valid for linear systems, while the system we consider here is nonlinear, we believe that this approach will give a good first order approximation especially if we consider small changes in the recharge.
To be able to use convolution for estimating the lens response to variable recharge, an impulse response has to be determined.This impulse response function defines the reaction of the system to a unit impulse, i.e. an impulse for which the integral over time equals one.Instead of using an impulse input, we used a step function.Since the impulse is the derivative of the step function in time and the system is considered to be linear, the impulse response function will be given by the time-derivative of the response to the step function.We have used a 0.5 mm day −1 change in the recharge to establish the step response functions.To minimize the effect of nonlinearities, both Introduction

Conclusions References
Tables Figures

Back Close
Full  3, the step response function s(t) can very well be fitted by an exponential function of the form The impulse response function is then given by the derivative of the step response function For a variable recharge P s in time, starting from an initial condition at time t = 0, namely the steady state lens with constant recharge P , the response of the lens in terms of the lens thickness in the middle of the field is then given by the convolution integral where the recharge P s needs to be given in mm day −1 and the factor 2 is a consequence of the determination of the step response for a change in recharge of 0.5 mm day −1 .
For a sinusoidal recharge pattern we can now obtain analytical expressions for the amplitude Q and the delay T of the sinusoidal response of the lens as is shown in the Appendix.Introduction

Conclusions References
Tables Figures

Back Close
Full

Volume change of the lens
We define lens volume as the volume of soil pores filled with fresh water (saturated and unsaturated) for a slab of 1 m thickness perpendicular to the 2-D-flow domain.To obtain a quantification of lens responses for thin lenses to recharge that varies as a function of time, we investigated relations between the most important parameter groups.To this extent, we simulated a large number of sinusoidal recharge variations (P s ) (Table 3).
For a broad range of parameters (Table 3), we found responses that nearly linearly relate the amplitude (A) and frequency (f ) to the volume variation of the lens.Volume variation is represented by a normalized volume deviation is the maximum and V the average lens volume.For a designated Rayleigh number (R), V n is linearly related to the period 1/f Ps where f Ps is the dimensionles group representing the frequency of P s (Eq.6a).
In Fig. 5a is shown how an increasing period leads to a larger volume variation.Differences in slope are primarily caused by different mass flux ratios (M), since variations are relatively larger for thinner lenses, which have a larger M. Volume variation is to a lesser extent increasing with larger recharge amplitudes represented by dimensionless group A Ps (Eq.6b).Multiplying f Ps and V n and relating this to M, results in Fig. 5b, where the slope now mainly depends on A Ps .The results of has a minor influence (less than 10 %), if we consider permeability ranges that seem plausible for drained agricultural fields.Equation ( 14) shows that the deviation of lens thickness from its average is linear and positively influenced by larger seepage and larger recharge amplitude and the average thickness itself.Higher frequency, on the other hand, diminishes the deviation of lens response.This is consistent with the findings of Cartwright (2004), which indicated that higher amplitude and lower frequency enhance each other in their positive effect on lens deviation.
Although Eq. ( 14) is an emprical relationship, the advantage is that it describes the impact of quite a number of parameters (namely 11) over a wide range.The normalized lens deviation does not increase linearly with increasing MA Ps for very large values, as is shown in Fig. 6b.For such conditions, saline water is being evaporated at the top boundary (or transpired by plants) and this decreases the response of the fresh water volume.The physical situation represented by these conditions implies presence of salt at the soil surface, and consequently significant amounts of salt in the root zone.
This will reduce primary production to a large extent.
The average lens volume V calculated for the numerical simulations for which Eq. ( 14) holds, deviates less than 5 % from the steady state lens that is formed when P s is replaced by a constant recharge equal to P s .As we have shown earlier (Eeman et al., 2011), the center of the mixing zone of steady state lenses is well approximated by the analytical model proposed by Maas (2007).Figure 7 illustrates that the relation between the total lens volume and its thickness in the middle of the field is nearly equal for the numerical simulations and the approximation by Maas.Therefore, we can replace V in Eq. ( 14) by V M , which can be calculated using the approximation of Maas ( 2007 in combination with his assumption of an elliptic lens shape this leads to a (half) lens volume (remember: our model domain is only half a lens in view of symmetry) Combining Eqs. ( 16) and ( 15), we obtain for the lens response based on recharge data and field characteristics For lenses with an average thickness of more than 3 m the approximation of Maas overestimates the lens volume (Fig. 7) derived from its thickness.This can be largely attributed to the outflow face in the model used by Maas, which creates a wider lens near the outflow region compared to the ditch outflow we model.For the lenses of interest here, as elaborated in Sect.2.3, the difference between the volume-thickness relations of the numerically and analytically calculated steady state lenses is less than 5 %.Therefore we can determine either and transform as required.

Thickness of the mixing zone
Since the mixing zones appear thick compared to total thickness of the fresh water lens when recharge varies as a function of time (Fig. 1b, De Louw, 2011), an estimate of its average position is not sufficient.Temporay saline water in the root zone may be caused by a thick transition zone, even when the average lens thickness covers the root zone.
We propose an analogy to a mixing zone that forms during the uniform motion of an initially sharp front through a porous medium.The width of the mixing zone for such a front depends on the distance it has traveled since t = 0.For a Fickian dispersion/ mixing regime, the standard deviation of a normally distributed concentration change is related to the diffusion/dispersion coefficient and time according to where a sharp interface at t = 0 is assumed.The standard deviation σ of the center of the mixing zone ( z), is calculated from the central moment of the concentration change distribution (inset Fig. 1a).For a lens, the assumption of an initially sharp front is not met; however, we can derive a relation between the distance traveled in a certain period and the width of the mixing zone.First, we define χ s as the traveled distance per sine-period of the center of the mixing zone ( z) in the middle of the field and |ν| is the average absolute velocity of this interface for sinusoidal recharge variation (Cirkel et al., 2012).Therefore and where A z is the amplitude of z, as Fig. 8 shows.We approximate hydrodynamic dispersion D by the longitudinal dispersivity multiplied by the average absolute interface velocity.This is justified since mixing in vertically expanding and shrinking lenses is dominated by longitudinal dispersion (Eeman et al., 2011) The combination of Eqs.(18-20) leads to a measure for the average variance of a mixing zone when recharge variation is a sinusoid, denoted as σ s 2 Because of the dispersion, there is no time t = 0 with a sharp interface as assumed for Eq. ( 18).Therefore, mixing zone width will partly depend on its initial thickness and position.Simulations show, however, that the relation between σ s 2 and 2α l χ s f is quite linear, as is illustrated in Fig. 9 for different mass flux ratios (M) and recharge curves (P s ).An explained variance for Eq. ( 21) of 0.85 is found for the smallest simulated mass Introduction

Conclusions References
Tables Figures

Back Close
Full The different coefficients for the different mass flux ratios (M) can be explained by the increasing total flux in these simulations.M = ρ max S/ρ 0 P s (Table 1) is varied by changing seepage S, which does not influence any of the other dimensionless groups.However, the total flux in the flow domain increases.This would lead to a wider mixing zone because larger fluxes lead to larger velocities.This mechanism is suppressed by streamlines that converge closer to the ditch where also velocity increases, as described in Eeman et al. (2011).
The approximation of net rainfall with a sinus ignores all irregular behavior (as is apparent in Fig. 3), and its suitability for this purpose therefore needs to be shown.Hence, we determined the thickness of the mixing zone also for a 15 yr period of daily recharge data in De Bilt (http://www.knmi.nl/datacentrum)compared to a sinusoidal recharge curve P s with the same average recharge, Fig. 10a.
The traveled distance (χ n ) of the center of the mixing zone ( z) of the simulation using daily weather data is calculated by The average variance of the lens when daily (i ) recharge records are used, σ n 2 , is calculated by Assuming that the variance depends on the traveled distance for both the sinusoidal and the natural recharge signals, using annual frequencies (f ) for both signals, we combine Eqs. or In Fig. 10a, we show net recharge patterns that were used as input.Figure 10c shows the center of the mixing zone z and Fig. 10d the associated mixing zone variance σ 2 .In Fig. 10b the dots indicate the ratio according to Eq. (24b) for individual years, for which significant deviations from the expected ratio of 1 are calculated.This can be attributed to erratic weather in terms of relatively dry or wet years and how these affect the initial conditions for the next year (Fig. 10c and d).The average ratio over a period of 15 yr is 1.01 to 1.03, which is a good indication that on average, the traveled distance is proportional to the mixing zone thickness.This method is therefore suitable to estimate the mixing zone thickness for any rainfall pattern or for example the effects of irrigation.The only condition is that a sufficiently long period is used to minimize the effect of initial conditions.To estimate the width rather than the ratio obtained by Eq. ( 24), a numerical simulation would be needed to establish a reference situation, from which all variations can then be derived analytically.

Delay and amplitude of lens response
The impulse response function was derived to apply convolution for thin lenses.We first derived the step response functions for a range of M and R, for small changes in constant recharge (0.5 mm day −1 ), as shown for a few example curves in Fig. 11, where the total gain ∆z and the shape parameter a were fitted according to Eq. ( 9).
This leads to Fig. 12a, b in which a and ∆z show a strong relation with M, and R has a limited influence on ∆z and a negligible impact on a.Note that Fig. 12 is not limited to sinusoidal recharge curves: it only relates a change in lens volume created by a change in constant recharge to differences in M and R.
For sinusoidal recharge curves, we can relate these parameters to the amplitude of the lens and the delay of the lens response compared to the recharge sinus, using the

HESSD Introduction Conclusions References
Tables Figures

Back Close
Full derivation as found in the appendix.The delay of a lens, calculated from Eq. ( 13b) was found to be always approximately 25 % (±2 %) of the period of the recharge variation (1/f ), irrespective of M,R,A Ps and f Ps .This is expected and understandable.The lens approaches its maximum when the recharge curve approaches the point where it changes from positive to negative.The maximum is reached slightly before this point, when recharge equals outflow.The little recharge still occurring is required to "maintain" lens thickness, due to the outflow into the ditch.The same is found when recharge changes from negative to positive: the minimum lens thickness is reached shortly after recharge has changed from negative to positive values.We compared this to the delay of the simulations used to establish Eq. ( 14), which leads to the same value, although the spread is slightly larger: 25 % (±3 %).Delays were determined for at 1 to make sure Eq. ( 13) is valid (as elaborated in the appendix).
Lens amplitudes determined from Eq. (13a), scaled with the average lens thickness, were compared to the numerical simulation results analysed in Sect.3.1, again for at 1. Differences are less than 5 %.Lens amplitude is again only slightly influenced by R. The relations with M,A Ps and f Ps are shown in Fig. 13.A thinner lens (larger M), is more influenced by changes in A and f .The enhancing effects of larger A and smaller f (Cartwright et al., 2004) is again confirmed.Only thin lenses are affected to an extent that may cause root zone salinization, with thickness changes of more than 30 % for not extreme values of A and f (Fig. 13c).Very clear from this analysis is that higher frequencies can only influence lens thickness significantly when amplitudes are extreme.

Conclusions
In this paper, we addressed the impact of temporal variations of net recharge at the soil surface on the thickness of, respectively fresh water lenses and their mixing zone at the fresh/salt interface.We investigated the impact with numerical 2-D simulations, varying dimensionless groups that follow from the basic governing flow and transport Introduction

Conclusions References
Tables Figures

Back Close
Full equations.The changes in both fresh water lens thickness and volume, and the mixing zone thickness, were related with simple linear functions of these dimensionless groups.
An empirical relation that was developed concerns the volume deviation of a thin fresh water lens from its average, in response to sinusoidal recharge variations (Fig. 7 and Eq.14).This relation clearly shows the positive effect of recharge amplitude and the negative effect of recharge frequency on lens thickness variation.Because the average lens thickness is hardly influenced by the amplitude and frequency of the recharge variation, the relation can be combined with a steady state approximation of lens volume under saline seepage conditions (Maas, 2007), Eq. ( 16).This relation holds for parameter combinations appropriate for a combination of realistic parameter values, and provides a simple, computationally fast estimate of the maximum and minimum lens thickness (Eq.16).
A first order approximation of the impulse response function for a thin lens was derived (Eq. 10), for which the parameters can be obtained from Fig. 12.With the convolution integral (Eq.11) it is then simple to determine a first order approximation of the position of the mixing zone for arbitrary recharge variation in time.For sinusoidal recharge patterns, amplitude of the lens with respect to the recharge changes can be easily calculated (Eq.13), whereas the delay turns out to be approximately 25 % of the sinus period, independent of lens conditions.Both the shallowest position and the time of occurrence can be determined.The former indicates the rooting depth at which plants may take up saline water, the latter provides the moment during the growing season that saline water reaches this minimum depth.This combination can be used to assess possible crop damage.Results obtained by convolution results are in very good agreement with numerically simulated results, which indicates that the approximation is very useful, even though the considered system is non-linear.
The average mixing zone thickness over a longer period (>10 periodic recharge cycles) can be estimated from the traveled distance of the average mixing zone position, (Eqs.21-24), for any recharge pattern.At least one reference simulation is required For larger values of t, i.e. if at 1, the terms with e −at will become negligible, leading to:

Conclusions References
Tables Figures

Back Close
Full Defining ε = arctan 2πf a which gives: we arrive at the following expression for the lens response to sinusoidal recharge: where are respectively the amplitude and the delay of the lens response to a sinusoidal recharge pattern.Introduction

Conclusions References
Tables Figures

Back Close
Full  Where is the average of , and the groups represent the dimensionless frequency and amplitude, respectively.In our analysis we will focus on these two groups, and the mass flux ratio and the Rayleigh number .Full  Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

aξµ
Shape parameter for step-and impulse-response function (−) f Frequency of recharge sinus (T −1 ) g Acceleration of gravity (L T −2 ) m sb Total salt mass flux across the ditch boundary (M L average abs.velocity of the center of the mixing zone in middle of field (L T −1 ) of the vertical center of salt mass change (L) A Amplitude of sinusoidal recharge curve (L T −1 ) A z Amplitude of the center of the mixing zone in the middle of the field (L) D m Apparent molecular diffusion coefficient, including the effect of tortuosity (L 2 T −1 ) D Hydrodynamic dispersion tensor, including apparent molecular diffusion (L 2 T −1 ) D eff Dispersion tensor, including D, D m , and the exchange between mobile/immobile water regions (V m , V resp.normalized, maximum and average volumes of fresh water (L 2 ) Z Maximum thickness of the fresh water lens in the middle of the field (L) α l , α t , α eff Longitudinal and transversal dispersivity (L) Conductance of the boundary between porous medium and ditch (L) φ m Fraction of water in the mobile phase (−) γ Constant in the equation of state for the fluid density (−) κ Intrinsic permeability of the porous medium (L 2 ) Fluid viscosity (M L −1 T −1 ) ρ, ρ 0 , ρ max Fluid density, fresh water density and maximum salt water density (M L −3 ) ∆z Ultimate gain for step and impulse response function (L) ∆ρ = ρ − ρ 0 Fresh water-salt water density difference (M L −3 ) σ, σ S , σ n standard deviation of the center of the mixing zone in the middle of the field, resp.for a sinusoidal recharge function and actual weather data (L) τ Integration variable in the convolution integral (T) χ , χ S , χ n Traveled distance per period of the center of the mixing zone in the middle of the field, resp.for a sinusoidal recharge function and actual weather data (L T −1 ) ω, ω max Salt mass fraction and maximum salt mass fraction (−) 3 Theory Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | the first term on the right side presents the hydrodynamic pore scale dispersion in the mobile phase and the second term gives the extra dispersion caused by the exchange between the mobile and immobile phases.The pore water velocity has Discussion Paper | Discussion Paper | Discussion Paper |a large influence, and is related with the mobile water fraction according to of this exchange process compared to the mobile phase dispersion is shown in Fig.5for r agg from 0.05 to 0.25 m (with D agg = 10 with r agg = 0.1 m).We show, for illustration, parameter combinations that lead to α eff = D eff /v = 0.1 m.Average recharge P is in the order of 0.3-3 mm day −1 .The combinations of parameters leading to α eff = 0.1 m give plausible values for aggregate size r agg , aggregate diffusivity D agg and mobile phase fraction φ m , whereas α eff = 0.25 leads to more extreme values.In view of all above, we use a longitudinal dispersivity of 0.1 m, and a transversal dispersivity of 0.01 m.
Discussion Paper | Discussion Paper | Discussion Paper | a negative and positive change in recharge has been applied and the average of the absolute values of the responses has been used to determine the step response function.Results show that for simulations with the M and R-values of Table Fig. 5 lead us to multiply A Ps and M to obtain one linear relation for the most important parameter groups, shown in Fig. 6a.The fitted linear relation is used to obtain an equation that relates lens deviation ∆V = V m − V to the average lens volume and three dimensionless groups β = 0.87 (95 % confidence interval from 0.86 to 0.88) is the slope of the line shown in Fig. 6a.The Rayleigh number R, which for practical situations mainly increases or decreases due to an increase resp.decrease in soil permeability κ, Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Table 1 :
Dimensionless groups.The inset figure shows the sensitivity of the steady state interface position for the groups considered by Eeman et al. (2011), with reference to a steady state lens of 5 m thick (their Fig. 6a).A positive change means a thicker lens.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 2 .Fig. 4 .Fig. 7 .Fig. 8 .Fig. 9 .
Fig. 1.(a) Representation of fresh water lens (white) floating on top of saline upward flowing groundwater (dark grey), with arrows that illustrate flow lines.The inset shows the concentration (solid) and concentration change (dashed) as a function of depth, where the latter indicates the position and width of the mixing zone.(b) Field measurements of resistivity with CVES (continous vertical electrical sounding) of lens and mixing zone, measured on Schouwen-Duiveland, The Netherlands, September 2007 (taken fromGoes et al., 2007; De Louw et al., 2011).The left part is a higher elevated sandy creek ridge, causing a thicker fresh water lens.The right part (110 to 140 m) is a thin lens on top of upward seeping saline groundwater, as schematized in (a).

Fig. 10 .Fig. 11 .Fig. 13 .
Fig. 10.(a) The recharge pattern of natural and sinusoidal recharge curves is shown.(b)The ratio between the traveled distance and variance (Eq.24b) for individual years (dots) and averaged over 15 yr (lines).(c) The center of the mixing zone for 15 yr of daily recharge (z 1 ) and (d) the associated variance (σ 2 n ) for lenses with a different average thickness.

Table 1 .
Eeman et al. (2011)s.The inset figure shows the sensitivity of the steady state interface position for the groups considered byEeman et al. (2011), with reference to a steady state lens of 5 m thick (their Fig.6a).A positive change means a thicker lens.

Table 2 .
Reference values of model parameters.

Table 3 .
Parameter variation of seepage (S), permeability (κ) and for the parameters of precipitation sinusoidal function (P s with P = 1, Eq. 5), for which the reference situation (bold values) is De Bilt (http://www.knmi.nl).