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

In coastal zones with saline groundwater, fresh groundwater lenses may form due to infiltration of rain water. The thickness of both the lens and the mixing zone, determines fresh water availability for plant growth. Due to recharge variation, the thickness of the lens and the mixing zone are not constant, which may adversely affect agricultural and natural vegetation if saline water reaches the root zone during the growing season. In this paper, we study the response of thin lenses and their mixing zone to variation of recharge. The recharge is varied using sinusoids with a range of amplitudes and frequencies. We vary lens characteristics by varying the Rayleigh number and Mass flux ratio of saline and fresh water, as these dominantly influence the thickness of thin lenses and their mixing zone. Numerical results show a linear relation between the normalised lens volume and the main lens and recharge characteristics, enabling an empirical approximation of the variation of lens thickness. Increase of the recharge amplitude causes increase and the increase of recharge frequency causes a decrease in the variation of lens thickness. The average lens thickness is not significantly influenced 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 centre 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 lowlying coastal zones, saline water is often present at a 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 S. Eeman et al.: Response to recharge variation of thin rainwater lenses reasons, which have caused subsidence of these low-lying areas and, therefore, increased the saline ground water pressure ( Van de Ven, 2003).For example, studies by Zeidler (1997), Day et al. (1995), Habibullah et al. (1999) and Osterkamp et al. (2001) conclude that salinisation in low-lying 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 minimise adverse effects of salinisation, water management and land use planning should be based on a proper understanding of the behaviour of thin fresh water bodies on saline groundwater and their response to temporal variability in recharge, whether seasonal, longer term, natural or human induced.
Related to the behaviour 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(ses) 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., 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 variability 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 at al. (1992), they found that the effect of the diurnal tidal pulse (0.5 to 2 m sea level rise) did not influence the 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. (2004) investigated much smaller spatial (horizontal and vertical) and temporal scales than Underwood et al. (1992) and Kiro et al. (2008), and neglected thickness variability 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 leads 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 variability.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 the 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 used to determine under which circumstances water quality is affected in the root zone.Lens properties that describe the depth at which saline water is found, are the lens thickness and the vertical extent of its mixing zone.We analyse the effect of recharge variations with different duration and intensity on lenses with a different thickness and mixing zone thickness, using numerical simulations.The first step is to establish a relation between recharge variation and lens thickness.Secondly, the effect of variation on the mixing zone is investigated, and finally delay and amplitude of lens systems with different recharge variations is calculated.Together these aspects provide a useful tool for a fast analysis of lens variability under a given recharge pattern.

Flow domain and system equations
The modelled 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 (Fig. 1a).
Recharge with fresh water (uncolored) occurs at the upper boundary, while vertical upward seepage of saline water (grey-coloured) 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.Lens thickness is defined as the distance between soil surface and the centre of the mixing zone, in the middle of the field (see also Sect. 2.3) 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 = xz).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 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 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 in the nomenclature).The mass flux ratio M is the ratio between the mass fluxes of the lower (seepage) and upper (net recharge) boundary, the Rayleigh number R is the ratio of the characteristic density induced flux and the average discharge P .
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) f P s = f Ln P and 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 Schematisation 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).SU-TRA (Voss and Provost, 2008) was used to carry out densitydependent numerical groundwater simulations, accounting for flow in the unsaturated zone.Spatial discretisation 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 discretisation 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 1.A longitudinal dispersivity of 0.1 m and a transversal dispersivity of 0.01 m were chosen.This choice is elaborated in Appendix A. Unsaturated soil parameters for a clay loam soil were used for the flow in the unsaturated zone.We tested the effect of different soil types on the lens thickness.This effect was negligible because the unsaturated zone is thin, and most of the time rather wet.For the initial conditions, steady state lenses were used that were obtained with a constant recharge rate of 1 mm d −1 .The volume and thickness of rainwater lenses and the thickness of the mixing zone for the numerical simulations were quantified using spatial moments (Eeman et al., 2011).Spatial moments efficiently summarise 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 modelled domain, the total volume of fresh water in the lens and transition zone can be inferred.The spatial moments of the vertical derivative 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 fluctuations, we varied the lens characteristics and the recharge.We varied the average lens thickness by varying M and R (Table 2), which leaves all other groups constant (using reference parameter values as given in Table 1).Only C, the dimensionless conductance of the soilditch interface, changes as we vary R. Since the influence of C is very small (inset Nomenclature), 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.We use the first spatial moment of dc/dz to determine lens thickness (see inset Fig. 1a).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 parameterise 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. 3 and the bold values in Table 2. 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 d −1 , while the period was varied from 1 yr to 1 week (Table 2).Although an annual frequency is most often expected, higher frequencies were simulated for    demonstrative purposes and because they have a resemblance with short-term precipitation events.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 weather events the mixing zone that we will propose (Sect.3.2).

System delay using convolution
To investigate the delay between the variations 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 variations 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 minimise the effect of nonlinearities, both a negative and positive change in recharge have been applied and the average of the absolute values of the responses has been used to determine the step response function.Results, as will be shown in Sect.4.3, show that for simulations with the M and R-values of Table 2, 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 Appendix B.

Volume variation 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 2D-flow domain.We use the zeroth moment (see Sect. 3.2) to calculate this volume.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 numerically simulated a large number of sinusoidal recharge variations (P s ) (Table 2).For a broad range of parameters (Table 1), 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 normalised volume deviation V n = (V m − V )/ V , where V m 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 P s where f P s is the dimensionless group representing the frequency of P s (Eq.6a).In Fig. 4a, it 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 P s (Eq.6b).Multiplying f P s and V n and relating this to M, results in Fig. 4b, where the slope now mainly depends on A P s .The results of Fig. 4 lead us to multiply A P s and M to obtain one linear relation for the most important parameter groups, shown in Fig. 5a.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 where β = 0.87 (95 % confidence interval from 0.86 to 0.88) is the slope of the line shown in Fig. 5a.The Rayleigh number R, which for practical situations mainly increases or decreases due to an increase resp.decrease in soil permeability κ, only has a minor influence (less than 10 %), if we consider permeability ranges that seem plausible for drained agricultural fields.Equation ( 12) 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 et al. (2004), which indicated that higher amplitude and lower frequency enhance each other in their positive effect on lens deviation.
Although Eq. ( 12) is an empirical relationship, the advantage is that it describes the impact of quite a number of parameters (namely 11) over a wide range.The normalised lens deviation does not increase linearly with increasing MA P s for very large values, as is shown in Fig. 5b.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, disappearance of an actual fresh water lens, 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. ( 12) 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 centre of the mixing zone of steady state lenses is well approximated by the analytical model proposed by Maas (2007).Figure 6 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. ( 12) 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. ( 14) and ( 13), 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. 6) 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.3.2, 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 et al., 2011), an estimate of its average position is not sufficient.Temporary 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 travelled 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 centre 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 travelled in a certain period and the width of the mixing zone.First, we define χ s as the travelled distance per sine-period of the centre of the mixing zone (z) in the middle of the field and |v| is the average absolute velocity of this interface for sinusoidal recharge variation (Cirkel et al., 2012).Therefore, where A z is the amplitude of z, as Fig. 7 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. ( 16)-( 18) leads to a measure for the temporal average variance (a suitable measure for thickness) 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. ( 16).Therefore, mixing zone width will partly depend on its initial thickness and position.Numerical simulations show, however, that the relation between σ s 2 and 2α l χ s f is quite linear, as is illustrated in Fig. 8 for different mass flux ratios (M) and recharge curves (P s ).An explained variance for Eq. ( 19) of 0.85 is found for the smallest simulated mass flux ratio, whereas for other mass flux ratios, the explained variance is larger than 0.96.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 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 behaviour (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. 9a.
The travelled distance (χ n ) of the centre of the mixing zone (z) of the numerical simulation using daily weather data is calculated by The average variance of the lens when daily (i) recharge records are used, σ s 2 , is calculated by Assuming that the variance depends on the travelled distance for both the sinusoidal and the natural recharge signals, using annual frequencies (f ) for both signals, we combine Eqs. ( 19)-( 21) to  In Fig. 9a, we show net recharge patterns that were used as input.Figure 9c shows the centre of the mixing zone z and Fig. 9d the associated mixing zone variance σ s 2 .In Fig. 9b the dots indicate the ratio according to Eq. ( 22b) for individual years, for which significant deviations from the expected ratio of 1 are simulated.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. 9c  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 travelled distance is proportional to the mixing zone thickness.This method of relating mixing zone thickness to the travelled distance of the recharge signal is, therefore, suitable for estimating 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 minimise the effect of initial conditions.To estimate the width rather than the ratio obtained by Eq. ( 22), 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 (derived from numerical simulations), for small changes in constant recharge (0.5 mm day −1 ), as shown for a few example curves in Fig. 10, where the total gain z and the shape parameter a were fitted according to Eq. ( 7).This leads to Fig. 11a and b in which a and show a strong relation with M, and R has a limited influence on z and a negligible impact on a.Note that     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 derivation as found in the appendix.The delay of a lens, calculated from Eq. (11b) was found to be always approximately 25 % (±2 %) of the period of the recharge variation (1/f ), irrespective of MRA P s and f P s .This is expected and understandable.The maximum lens volume is reached at the moment   when recharge equals outflow during a phase of declining recharge.The same is found for the minimum lens thickness, which is reached when recharge equals discharge during a phase of increasing recharge.We compared this to the delay of the simulations used to establish Eq. ( 12), which leads to the same value, although the spread is slightly larger: 25 % (±3 %).Delays were determined for at 1 to make sure Eq. ( 11) is valid (as elaborated in Appendix A).
Lens amplitudes determined from Eq. (11a), 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 MA P s and f P s are shown in Fig. 12.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 salinisation, with thickness variations of more than 30 % for not extreme values of A and f (Fig. 12c).It is clear from this analysis that higher frequencies can only influence lens thickness significantly when amplitudes are extreme; their main effect is on the thickness of the mixing zone, as explained in Sect.4.2.

Discussion and 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 equations.The variations 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. 5 and Eq.12).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. ( 13).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.15).
The average mixing zone thickness over a longer period (< 10 periodic recharge cycles) can be estimated from the travelled distance of the average mixing zone position, (Eqs.19-22), for any recharge pattern.At least one reference simulation is required to derive an estimate of the absolute value of the thickness for any sinusoidal variation of recharge.The mixing zone analysis shows that the influence of short-term precipitation events on the thickness of the mixing zone is significant, in spite of their limited effect on lens thickness.
A first order approximation of the impulse response function for a thin lens was derived (Eq.8), for which the parameters can be obtained from Fig. 11.With the convolution integral (Eq.9) 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 variations can be easily calculated (Eq.11), 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 are in very good agreement with numerically simulated results, which indicates that the approximation is very useful, even though the considered system is nonlinear.
Together the proposed equations provide a first order approximation for all aspects of interest concerning salinity for thin fresh water lenses on saline, upward seeping ground water.The attraction of our analysis is nicely illustrated by considering the impact of more realistic erratic rainfall (measured daily time series) on mixing zone thickness.Whereas most of the analysis was for regular sine variation, also for erratic rainfall, the mixing zone thickness and its variation at first order are well reproduced.The approach will be most successful when based on one or a few numerical calculations to establish a good set of reference parameters for the specific field and soil type (including parameters for the unsaturated zone), especially to determine the thickness of the mixing zone.The approach can be used to develop a tool for upscaling to, for example, a more regional analysis of salt sensitivity of agricultural soils.For such an analysis, additional geological and geographical parameters, in particular, the layering of soil (De Louw et al., 2011) and spatially different drainage levels may have to be accounted for.

Appendix A 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 et al., 1992;Kaleris, 2006;Kaleris et al., 2002).The value of the dispersivities that follows from Fiori and Dagan (2000) would easily be 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. A1, 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 A1a 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   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)    +a 2 and T = 1 2πf arctan 2πf a (B7 or 11a, b) are, respectively, the amplitude and the delay of the lens response to a sinusoidal recharge pattern.

Fig. 1 .
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 (continuous vertical electrical sounding) of lens and mixing zone, measured on Schouwen-Duiveland, The Netherlands, September 2007 (taken fromGoes et al., 2009; 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 schematised in (a).

Figure 3 :Fig. 2 .
Figure 3:Fig.2. Sketch of the geometry and boundary conditions for a half field and its drainage area (ditch).The black lines are the domain boundaries.Areas 1 to 4, separated by grey lines, represent different discretisation zones.Length and width are indicated and the ditch has a triangular cross-section and a water level of 1 m below soil surface.

Fig. 3 .
Fig. 3. Daily weather data from De Bilt 2008, The Netherlands (KNMI) and a recharge sinus (fitted with average weather date from 1971-2000) used for numerical simulations.

Fig. 4 .
Fig. 4. (a) Linear relation between dimensionless recharge period 1/f P s and the normalised volume deviation V n for P = 1 mm day −1 .(b) Multiplying f P s and V n on the y-axis and relating this to the mass flux ratio (M) shows a linear relation dependent on A P s .

Fig. 5 .
Fig. 5. (a) Relation between mass flux ratio M times A P s and V n f P s for R = 11.87.The red line shows the fitted linear model, y = 0.87 x+ 0.18 with an explained variance of 0.99.(b) Indication that for larger values of MA P s the relation is nonlinear.

Fig. 6 .
Fig.6.The relation between lens thickness and lens volume according to the assumption an elliptical lens shape(Maas, 2007, his   Eq.5a-c, compare this paper Eq.16) and the numerical simulations.

Fig. 7 .
Fig. 7. Position (z) of mean interface as a function of time (t) for sinusoidal change, in support of Eq. (19).

Fig. 9 .
Fig. 9. (a) The recharge pattern of natural and sinusoidal recharge curves is shown.(b) The ratio between the travelled distance and variance (Eq.22b) for individual years (dots) and averaged over 15 yr (lines).(c) The centre 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.

Fig. 10 .
Fig. 10.Change of lens thickness for a 0.5 mm −1 day change of a constant recharge, calculated for thin lenses with different M and R. Fitted with exponential Eq. (7). a determines the shape of the curvature and Z is the ultimate value of the thickness change.

Fig. 11
Fig. 11 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 derivation as found in the appendix.The delay of a lens, calculated from Eq. (11b) was found to be always approximately 25 % (±2 %) of the period of the recharge variation (1/f ), irrespective of MRA P s and f P s .This is expected and understandable.The maximum lens volume is reached at the moment 31 Figure 11:

Fig. 12 .
Fig. 12. Contours of change in lens thickness as a fraction of the average lens thickness as a function of frequency and amplitude.(a) M = 0.51, (b) M = 2.03, (c) M = 5.07.For the reference situation this leads to lenses with a thickness of 4.8, 2.3 and 1.1 m, respectively.

Fig. A1 .
Fig. A1.Profiles of salinity as a function of depth for thin lenses.Numerical simulations with different dispersivities are compared to several representative soil profiles measured on two field sites on the island of Schouwen-Duiveland, The Netherlands.
eff = Dφ m + (1 − φ m ) r 2 agg v 2 15 D agg (A1)where 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 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.A2for r agg from 0.05 to 0.25 m (with D agg = 10 −10 m 2 s −1 ) and D agg from 5 × 10 −11 to 5 × 10 −10 m 2 s −1 (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 the above, we use a longitudinal dispersivity of 0.1 m, and a transversal dispersivity of 0.01 m. 32 Figure AI2:

Fig. A2 .
Fig. A2.Relations based on Eq. (AII1) between the mobile fraction (φ m ), average net precipitation P and aggregate diffusivity D agg for a given aggregate radius (r agg ) of 10 −2 m (a) and with aggregate radius (r agg ) for a given aggregate diffusivity (D agg ) of 10 −9 m 2 s −1 (b), leading to an effective dispersivity (α eff ) of 0.1 m.