Spatial distribution of solute leaching with snowmelt and irrigation : measurements and simulations

Introduction Conclusions References Tables Figures


Introduction
Groundwater contamination by nutrients or chemicals will be enhanced by preferential flow in the unsaturated zone.Preferential flow can be caused by macropore flow (Jarvis, 2007), by small scale differences in hydraulic properties (Roth, 1995), or by water repellency (Van Dam et al., 1990), amongst others.To account for preferential flow in modelling heterogeneous soils, several approaches exist, which were reviewed by Feyen et al. (1998) and Šimunek et al. (2003).Roth (1995) included soil heterogeneity in a numerical model, by the use of a random scaling factor for the retention curve and the saturated hydraulic conductivity, where heterogeneous water flow was studied.Roth and Hammel (1996) extended the study by including solute transport.Roth (1995) showed that the infiltration rate determines which parts of the soil will transport most water and solutes.Another modelling approach is the use of independent stream tubes, which each have a different velocity and dispersion coefficient (Vanderborght et al., 2006;Russo and Fiori, 2009).Solute transport can evolve from a stochastic-convective (independent stream tubes) to a convective-dispersive regime with increasing depth (Seuntjens et al., 2001).
To include soil heterogeneity, a random distribution is needed for the spatial variability of the hydraulic properties like the saturated hydraulic conductivity.This distribution can be based on many samples of a soil profile, on which the hydraulic properties are determined (Rockhold et al., 1996;Hammel et al., 1999).Both used a scaling factor for similar media to define the random distribution of the hydraulic properties (Miller and Miller, 1956;Warrick et al., 1977).
Besides numerous measurements of the soil hydraulic properties, multi-compartment samplers (MCS) can be used

D. Schotanus et al.: Spatial distribution of leaching
to investigate the effect of heterogeneous flow on solute leaching of undisturbed soils in the field (Holder et al., 1991;Quisenberry et al., 1994;Bloem et al., 2009Bloem et al., , 2010)).The instrument has a porous plate to which suction is applied, and it consists of several cells in which drainage is collected.The volume of the drainage, and the concentration of a tracer or reactive solute are measured in each cell, to quantify the spatial variability in solute leaching.The step from this type of experiments to models that quantify solute leaching under different conditions, has not been made yet.
We use results from experiments with an MCS for transport modelling in a heterogeneous soil, by basing the extent of soil heterogeneity in the model on these experiments.To quantify the spatial variability in solute leaching, and to investigate the effect of infiltration rates on this spatial variability, we performed two field experiments with an MCS in exactly the same location (Schotanus et al., 2012).The first experiment was done during snowmelt, with high infiltration rates.The second experiment was done with irrigation, with lower infiltration rates.In the snowmelt experiment the spatial differences in the solute concentrations were larger than in the irrigation experiment.This is possibly due to lower lateral exchange with a shorter residence time, which results in less dilution, and larger differences in the concentrations.In the irrigation experiment, more isolated peaks in the concentration were found (spatial autocorrelation of 0.26 versus 0.61 during snowmelt), from which can be concluded that heterogeneous flow in the soil was caused by small differences in the soil hydraulic properties.We use these experiments to generate a random field for the scaling factors.This field is then used in a model to further investigate the effect of the flow rate on the leaching of a tracer and a degradable contaminant.Furthermore, the effect of snowmelt on solute leaching is studied.
A tool to study solute leaching simultaneously in space and time is the leaching surface (De Rooij and Stagnitti, 2002).With a leaching surface, the scaled solute flux density is plotted for each cell of a sampler or model, and the cells are sorted descendingly by total leached mass.Bloem et al. (2008) applied the leaching surface to flow in a heterogeneous aquifer.The leaching surface has not yet been used for degradable solutes.By comparing leaching surfaces of a tracer and of a degradable solute, the effect of a transient flow rate on the heterogeneous leaching of a degradable solute can be investigated.
As experiments with multi-compartment samplers are rarely modelled, and as the leaching surface was not applied yet to degradable solutes, our objectives are to: (1) develop and test a new approach, using multi-compartment sampler data, on which the random field for the scaling factor is based; (2) investigate the effect of the influence of the infiltration rate on the leaching of degradable solutes; (3) apply the leaching surface to degradable solutes.

Experiment
A multi-compartment sampler (MCS) (Bloem et al., 2010) was installed at the field station Moreppen, near Oslo Airport, Norway (French et al., 1994).The field station is located in a flat area with coarse glaciofluvial sediments (sand and gravel) (French and Van der Zee, 1999).The soil consists of 15 % fine sand, 75 % medium and coarse sand, and 10 % gravel (French et al., 1994).The saturated hydraulic conductivity is 6.65×10 −4 m s −1 (French et al., 2001).The soil surface was covered with short grass.The size of the MCS is 31.5×31.5 cm 2 , and consists of 100 separate drainage collectors.The surface of the MCS consists of porous metal plates, to which pressure can be applied.Technical details about the MCS can be found in Bloem et al. (2010).
From a trench a horizontal tunnel was dug, leaving the soil above this tunnel undisturbed.The MCS was installed at 68 cm from the trench wall, and at 51 cm below the soil surface (Fig. 1).A 2 mm thick layer of wetted soil from the tunnel was applied to the surface of the MCS, to ensure a good contact between the MCS and the soil above it.After installation, the tunnel was backfilled, to avoid boundary effects.Four tensiometers were installed near the trench wall, at the same depth as the MCS.The average pressure head, which was measured by these tensiometers was applied to the MCS, plus 15 cm extra pressure head to compensate for a pressure head drop in the porous metal plates (Bloem et al., 2009).The pressure in the MCS varied in time.
Two experiments were done, while the MCS remained in the same location: one during snowmelt (26 March-23 May 2010) and one with irrigation (23 May-4 July 2010).For the snowmelt experiment 1092 g m −2 propylene glycol (PG) and 10 g m −2 bromide was diluted in 2 L m −2 , and sprayed homogeneously on top of an undisturbed snow cover (26 March 2010).De-icing fluid (Kilfrost, 2012) containing PG was diluted to reach this applied mass.The application area was sufficiently large, such that boundary effects can be ignored.No ice layer was observed under the snowcover, on top of the soil surface.Thus, no ponding could occur on top of the soil surface.
For the irrigation experiment 1103 g m −2 PG and 10 g m −2 bromide was diluted in 5 L m −2 , and sprayed homogeneously on the soil surface at 23 May 2010.In between the irrigations, the soil surface was covered with plastic, to prevent evapotranspiration and the infiltration of rainwater.By means of a nozzle, the soil surface was irrigated with tapwater, aside from 12 June.This was a rainy day and 7.8 mm of rainwater infiltrated, while the plastic was removed from the soil surface.On 13 June, 7.5 mm water was irrigated.The irrigation rate was around 6.6 mm h −1 .This was low enough, such that neither surface runoff, nor ponding was observed.Between 31 May and 2 June, evaporation was measured in a pan under the plastic.The measured evaporation was 1 mm d −1 .These were warm days, the average evapotranspiration rate during the entire irrigation experiment probably was lower.Drainage was stored in the MCS, the sampling scheme depended on the amount of drainage.For the snowmelt experiment, samples were taken every day from 1 April until 15 April.Thereafter, samples were taken at 17 and 21 April, and 9 and 23 May.During the irrigation experiment, samples were taken every second day, on the same day as, and prior to, the irrigations.After collecting, the volumes of the samples were measured by weighing.If a sample was smaller than 4 mL, it was too small to measure the bromide concentration, and it was stored cool.In the next sampling round, the new sample was added to this stored small sample.In samples larger than 4 mL, the bromide concentration was measured with an Orion 9635BNWP Bromide Ion Selective Electrode.Moreover, in the irrigation experiment, the PG concentration was measured with gaschromatography (Trace GC Ultra, Thermo Scientific).We assumed that degradation of PG during the snowmelt experiment was low, as the soil temperature was low, and the PG concentration in the applied solution was high (Jaesche et al., 2006).Therefore, the PG concentration was not measured in the snowmelt samples.
More details about the experiments can be found in Schotanus et al. (2012).

Models
To simulate the results of the experiment, we used Hydrus-2D, which is a model that can simulate water and solute transport in variably saturated porous media ( Šimunek and Sejna, 1999).We did not intend to simulate the results of the experiment such that the simulated and the experimental leaching were similar in exactly the same location in the horizontal plane.Instead, we were looking for a distribution of the soil hydraulic properties in a vertical plane, which will give the same spatial variability of drainage, throughout the entire experiment.The simulated drainage was sorted and cumulated, and then compared to the measured sorted and cumulated drainage.The degree of correspondence was determined with the mean absolute error.We used a model because the breakthrough curves (BTCs) of the experiment were truncated, because the time for the experiment was limited.From the modelling we can use the complete BTCs, and thus study the leaching of a degradable solute and a tracer in more detail.
For the saturated hydraulic conductivity and the pressure head, the standard deviation, and the correlation lengths in depth and width direction were determined with Hydrus-2D, using the spatial distribution of the drainage.The parameter set resulting in the smallest deviation between observed and calculated drainage was selected.As a criterion for the selection of the most appropriate random distribution for this soil, the total drainage per cell was used, sorted in a decreasing order, and then cumulated.Random fields with Millersimilarity were generated in Hydrus-2D.For details on the generation of these fields, we refer to the Hydrus manual ( Šimunek and Sejna, 1999).A vertical plane of 2 m width and 1.5 m deep was used.The groundwater level was fixed at 1.5 m depth.The actual groundwater table is situated at 4 m depth.As the soil has a coarse texture, the groundwater did not influence the leaching at 0.5 and 1 m deep, which were the depths of interest in this study.Density driven flow was ignored.The degradation of propylene glycol depends on temperature (Jaesche et al., 2006), but this was ignored.Degradation was modelled as a first-order process, with a half-life time τ of 10 d.In the field, the half-life time of PG can be up to 17 d (French et al., 2001), however, the half-life  Fig. 2. Precipitation or infiltration rate (a).1997 with snowmelt, solute application with a pulse of 6 days, snowmelt from day 1 to 6, length of weather series is 367 days (b).1997 without snowmelt, solute application with a pulse of 1 day, length of weather series is 244 days.The lines indicate equal parts of the precipitations series in (a) and (b).(c).snowmelt experiment 2010, solute application with a pulse of 1 day, length of weather series is 62 days.The day is the number of days since the first solute application.time can vary in time, for instance as soil temperature or wetness vary.Adsorption was not considered to occur in the model (French et al., 2001).
Output of concentrations and water flow velocities were given four times a day.In the comparison between the leaching surfaces of the snowmelt experiment and the simulation, the simulated values are aggregates to once a day, because the measurements were done daily.
The upper boundary was an atmospheric flux.To investigate the effect of the atmospheric flux on solute leaching, three different atmospheric input fluxes were used, measured at the weather station at the airport where the field station is situated: the snowmelt of the year 2010, and the year 1997, with and without snowmelt (Fig. 2).The year 1997 is a generally dry year, but had a thicker snowcover than the year 2010.For the year 2010, the infiltration rate was calculated from the snowmelt and precipitation, which was measured at the field site (Schotanus et al., 2012).For the year 1997, we made an approximation for snow formation and melt in a pre-processing routine: (1) where S is the snowdepth (m), b is a fitting parameter (m snow/m water), P is the precipitation (m), T is the temperature ( • C), M is the snowmelt, c is a fitting parameter (m • C −1 ), T melt is the critical temperature for snowmelt ( • C), B is the total storage available for water in the snowcover (i.e., porosity of the snowcover × depth of the snowcover, m), D is the depth of the snowcover (m).Fresh snow S is added to the snowcover D. When the temperature is higher than 0 • C, P and M are added to the fluid water depth W .When W exceeds B, this is called infiltration.The resulting infiltration rate was used as precipitation in Hydrus-2D.The parameters b and c were fitted with data from other years, and were 8 m snow/m water, and 0.001 m • C −1 , respectively.For the year 1997 without snowmelt, the high infiltration rate from the snowmelt, and the low infiltration rate during winter time were removed, and only the precipitation was used.
For the year 1997 with snowmelt, solute application started at the beginning of the snowmelt.Just before snowmelt, the soil is dry, as there is no infiltration during winter, because all precipitation is stored in the snowcover.For the year 1997 without snowmelt, there is infiltration throughout the year, as precipitation is not stored in the snowcover.Therefore, for the year 1997 without snowmelt, the soil moisture content is higher at the moment of solute application than for the year 1997 with snowmelt.In Fig. 2b the same infiltration rate is used as in Fig. 2a, but the day of solute application differs and, thus, the infiltration rate is shifted in the Figure .By using time periods with different colours, the shifting in time is made clear.

Data analysis
Moment's analysis is used to characterise the average transport of the solute plume (Burr et al., 1994).The first temporal moment of the concentration is the mean breakthrough time (Govindaraju and Das, 2007): where M 1 (z) is the mean breakthrough time (d), z is depth (L), which is 0 at soil surface, t is the time (d), and C is the solute flux concentration (M m −3 ).The leached mass is calculated as the sum of the convective and the dispersive flux (Kreft and Zuber, 1978): where LM is the leached mass (M), v is the water flux (m d −1 ), and D is the dispersion coefficient (m 2 d −1 ).
From the leached mass, the solute flux density is calculated (M L −2 T −1 ), which is the leached mass per area at a certain depth per time interval.To visualise the data from the experiments and the simulations, we use leaching surfaces (De Rooij and Stagnitti, 2002).In a leaching surface the leaching in space and time is shown simultaneously.The cells are sorted descendingly by the cumulated leaching per cell.This means that in the leaching surface the high leaching cells can be found at a low cumulative area, whereas the low leaching cells can be found at a high cumulative area.The solute flux density is scaled with the total leaching of all cells (M L −2 for the experiment, and M L −1 for the simulations).We use the total tracer leaching to scale the solute flux density of the degradable solute.In the snowmelt experiment, the bromide recovery was 43 %.In the irrigation experiment, the bromide recovery was 42 %, and the PG recovery was 32 % (Schotanus et al., 2012).De Rooij and Stagnitti (2002) use a 3-D plot for the leaching surfaces.We use 2-D plots with a colour scale, because 2-D plots are easier to interpret, and to compare with each other.
The leaching surface has two marginals: one in the time axis and one in the spatial axis.The marginal in the time axis is the breakthrough curve of the total area, and the marginal in the spatial axis is the spatial solute distribution.The spatial distributions of the solute leaching under different conditions can be compared by comparing the marginal in the x-axis of the leaching surfaces.The leaching surfaces will be compared with the following characteristics: the instant that the leaching starts (i.e., the amount of cumulative drainage until the first detectable concentration), the magnitude of the solute flux density, the tailing, and the fraction of the soil that contributes to solute leaching.In a cell or node that does not contribute to leaching, no drainage occurs, or the concentration is 0, both throughout the experiment.

Leaching surfaces from experiments
Figure 3 shows the 2-D leaching surface of bromide during the snowmelt experiment.To facilitate the comparison between the different experiments, the cumulative drainage since solute application is used as a time axis, instead of the number of days, in accordance with Wierenga (1977).In a leaching surface, the highest leaching cells can be found at a low cumulative area.The amount of leaching decreases with an increasing cumulative area.For the snowmelt experiment, the instant that the leaching in a cell starts, is generally later with a decreasing amount of leaching.Generally, the leaching per day decreased with a decreasing leached amount per cell.The experiment was stopped after 111 mm of drainage, which explains the lack of bromide leaching thereafter.About 85 % of the cells contribute to the leaching.The CV in the drainage of the cells was 0.9, the CV of the leached bromide was 1.1.
Figure 4a shows the leaching surface of bromide during the irrigation experiment.Here, the leaching of the high and mean leaching cells starts at the same time, after 50 mm drainage.The leaching starts after a larger amount of drainage than in the snowmelt experiment (20 mm), because the soil was wetter at the beginning of the irrigation experiment than of the snowmelt experiment (pressure head was −27 and −35 cm, respectively).The highest solute flux density observed in any of the cells is lower than in the snowmelt experiment (0.005 against 0.01 d −1 ), probably caused by the lower water flux in the irrigation experiment (6 mm d −1 during irrigation, and 16 mm d −1 during snowmelt on average).Fewer cells than in the snowmelt experiment contribute to the leaching, about 70 %.This is supported by the CV in the drainage of the cells, which was 1.1, higher than in the snowmelt experiment.The CV of the leached bromide was 1.3, also higher than in the snowmelt experiment.
Figure 4b shows the leaching surface of propylene glycol (PG) during the irrigation experiment.In most cells, the leaching of PG starts earlier than of bromide.After 50 mm of drainage, the solute flux density of PG is higher than of bromide, which is possibly caused by density driven flow.The density of pure de-icing fluid is 1.043 times the density of pure water (Kilfrost, 2012).After dilution the density of the applied solution was approximately 1.005 times the density of water.The tailing is less, probably due to degradation of PG.Slightly fewer cells contribute to the leaching of PG than of bromide, 70 %, which also can be caused by degradation.
The solute leaching becomes more homogeneous with an increasing infiltration rate (comparing Figs. 3 and 4a).If the heterogeneous leaching would be caused by macropores, the leaching would become more heterogeneous with an increasing infiltration rate (Jarvis, 2007).Therefore, it is concluded that the heterogeneous leaching is caused by small scale differences in the soil hydraulic properties.The soil heterogeneity can thus be described with Miller-similarity.

Model parameterisation
The observations from the irrigation experiment will be used to fit the parameters of the model, because in this experiment the atmospheric boundary condition is well known.Using the inverse mode in Hydrus-2D, the parameters α and n from the van Genuchten equations (van Genuchten, 1980) were fitted to the measured average pressure heads at 0.51 m depth, where the MCS was located.shows the observed and simulated pressure heads at 0.5 m for the irrigation experiment.The optimal value for α was 14.85 m −1 (95 % confidence interval: 14.38-15.30)and for n 3.165 (-) (95 % confidence interval: 3.09-3.24).The linear regression coefficient between the observed and fitted values (R 2 ) was 0.83.The values of the other soil hydraulic parameters were the residual water content θ r = 0.045 m 3 m −3 (from the category "sand", Carsel and Parrish (1988)), the saturated water content θ sat = 0.33 m 3 m −3 (measured), and the saturated hydraulic conductivity K sat =6.65 • 10 −4 m s −1 (French et al., 2001).
After fitting the soil physical parameters, the standard deviation and the correlation lengths in width x and depth z direction were varied manually for K sat and the pressure head h.The parameter set resulting in the smallest deviation between observed and calculated drainage was selected.As a criterion for the selection of the most appropriate random distribution for this soil, the total drainage per cell was used, sorted in a decreasing order, and then cumu-lated.Random fields with Miller-similarity were generated in Hydrus-2D.First, the spatial discretisation needed to capture all small scale processes was chosen.A random field with standard deviation σ = 0.5, correlation length in x-direction λ x = 0.05 m, and correlation length in z-direction λ z = 0.6 m was used to simulate the irrigation experiment, with spatial discretisations of 0.0125 and 0.025 m.The values for σ , λ x , and λ z were chosen such that they are extreme values, with a high standard deviation, and correlation lengths that lead to narrow and long flow channels.These values are only used to test the discretisation, for other simulations other parameter values will be used.The sorted cumulated drainage was similar for the discretisations of 0.0125 and 0.025 m.Therefore, we use a spatial discretisation of 0.025 m, as it captures the small scale processes, and saves computation time.Furthermore, it is concluded that a spatial discretisation of 0.025 m is small enough, when a correlation length of 0.05 m is used.This result is in contrast with Ababou et al. (1989), who found that the spatial discretisation should be at least four times smaller than the correlation length.As the values for σ , λ x , and λ z were extreme values, the conclusion that a discretisation of 0.025 m is sufficient, should also hold for smaller σ , larger λ x , and smaller λ z .The cell size of the multi-compartment sampler is 0.0315×0.0315m 2 , similar as the discretisation.
In a preliminary investigation of the most suitable parameters for these random fields, σ was 0.1, 0.25, 0.5 or log(0.15),λ x was 0.05, 0.1, 0.15, 0.2, 0.3 or 0.4 m, and λ z was 0.1, 0.15, 0.2, 0.3, 0.4 or 0.6 m.With these parameter values, random fields for the scaling factor were generated, which were used to simulate the irrigation experiment.Figure 6 illustrates the criterion for the observed sorted drainage and the results from two simulations.In this figure the mean absolute error between the observed and simulated cumulated and sorted drainage is 0.022 and 0.19 (-).
After the preliminary investigation, five parameter sets were selected (σ = 0.5, λ x = 0.05 m, λ z = 0.6 m; σ = 0.25, λ x = 0.05 m, λ z = 0.6 m; σ = 0.5, λ x = 0.3 m, λ z = 0.6 m; σ = 0.5, λ x = 0.05 m, λ z = 0.15 m; σ = 0.5, λ x = 0.1 m, λ z = 0.15 m).These were used to generate five random fields for each parameter set, to investigate the effect of a particular field, within a parameter set, on the selection criterion.The resulting twenty-five random fields were used to simulate the irrigation experiment.The mean absolute error between the observed and simulated cumulated and sorted drainage was calculated (Table 1).The mean absolute error depends on both the particular random field, and on the parameters of the random field.From Table 1 is concluded that σ = 0.5, and λ x = 0.05 m give results that correspond best with the observed data.The λ z does not influence the mean absolute error much.The fields with the two best parameter sets (σ = 0.5, λ x = 0.05 m, λ z = 0.6 m; σ = 0.5, λ x = 0.05 m, λ z = 0.15 m) were used to simulate the snowmelt experiment.Again, the mean absolute error between the observed and simulated cumulated and sorted drainage was calculated.Table 1 shows that σ = 0.5, λ x = 0.05 m, λ z = 0.15 m gives the smallest mean absolute error for the snowmelt experiment.
To further examine whether the parameterisation of the model is a good representation of the field site, first we compare the leaching surfaces of the measurements of the snowmelt experiment (Fig. 3), and the simulation of the snowmelt experiment (weather series of Fig. 2c). Figure 7 shows the leaching surface from the simulation of the snowmelt experiment with the realisation with the parameters σ = 0.5, λ x = 0.05 m, λ z = 0.15 m that had the smallest mean absolute error.Leaching starts after 15 mm of drainage.After 40 mm of drainage the highest solute flux density occurs in both the experiment and the simulation.To compare the magnitudes of the solute flux densities, the solute flux erated random fields with Miller-similarity.density of the simulation should be corrected with 59/100, as the simulation consists of 59 cells, and the experiment of 100 cells.The maximum solute flux density is slightly larger for the simulation (0.010 and 0.012 d −1 ).In the experiment 90 % of the soil contributes to solute leaching, against 98 % in the simulation.The mean absolute error between the observed and simulated cumulated and sorted solute leaching was 0.018 (-), which is smaller than the mean absolute error for the drainage for this random field (0.022).
Comparing Figs. 3 and 7, we conclude that the model captures the magnitude and moment of leaching of the experiment sufficiently well.Comparing the marginal of the spatial axis, which is the spatial distribution of the solute leaching, the model corresponds well with the experiment.Therefore, we conclude that the parameters σ = 0.5, λ x = 0.05 m, λ z = 0.15 m can be used to quantify the heterogeneity of this soil.The realisation with these parameters that had the smallest mean absolute error was selected for further simulations, to study the effect of different infiltration rates on the heterogeneous solute leaching in more detail.Figure 8 shows the random field for the scaling factor of the saturated hydraulic conductivity with σ = 0.5, λ x = 0.05 m, λ z = 0.15 m which is used in the simulations.
When the total leached amount of the cells in Fig. 7 is related to the scaling factor in the particular cells at 0.5 m depth (Fig. 8), the leached amount generally decreases with an increasing scaling factor.The saturated hydraulic conductivity (French et al., 2001) is high compared to the water flux (respectively 57 m d −1 versus 0.016 m d −1 ).As a result, the parts of the soil where the scaling factor is lowest, transport most water and solutes, as was shown by Roth (1995).Where the scaling factor is small, the saturated hydraulic conductivity is still high, because of the high mean saturated hydraulic conductivity, therefore, transport is fast in this soil.

Transient simulations
To investigate the effect of snowmelt on the leaching of a tracer and a degradable solute, a simulation with the weather series of the year 1997 with snowmelt (Fig. 2a) was done.
Figure 9a shows the simulated leaching surface for a tracer at 0.5 m depth with the random field of Fig. 8.The depth of the snowcover was 240 mm.In the highest leaching cells, the leaching starts immediately with the drainage.The magnitude of the solute flux density is higher as for the experiment and Fig. 7.The highest solute flux density occurs after 40 mm of drainage, similar as in the experiment and in Fig. 7.The highest solute flux density continues longer, because the snowcover is thicker, and the solute was applied during the snowmelt period, which was 6 days.The entire soil contributes to solute leaching during the snowmelt.After 210 mm, when all snow had infiltrated, there is some drainage with little solute.The solute still leaches, but due to the low drainage, the solute flux density is very small.Then, after a precipitation event, the water and solute fluxes increase again.In the first precipitation events since snowmelt (after 220 mm of drainage), the leaching occurs in 80 % of the cells.The percentage of the soil that contributes to solute leaching decreases with increasing time since snowmelt.In Fig. 9a the solute flux density is relatively high in the entire area, which larger than in Figs. 3, and 7.This can be caused by the higher infiltration rate in the year 1997, as the snowcover was thicker in 1997 than in 2010, this results in a higher water content of the soil.As a result, a larger fraction of the soil is highly conductive than in 2010.
Figure 9b shows the simulated leaching surface for a degradable solute with half-life time τ = 10 d.Until 210 mm (during snowmelt), the leaching of the degradable solute is similar as of the tracer, both in space and time.Also the magnitude of the solute flux density is similar.This means that the infiltration rate is high compared to the degradation rate of the solute.After 210 mm, the leaching of the degradable solute differs from the leaching of the tracer.After a precipitation event, at 220 mm drainage, the tracer leaches in approximately 80 % of the cells, while the degradable solute leaches in only 1 % of the cells.The reason for this is that the degradable solute is mostly degraded in the period between the snowmelt and the precipitation event, which is 32 days long.As a result, then the solute flux density of the degradable solute is smaller than for the tracer.
From Figs. 3, 7, and 9 can be concluded that the leaching surface is highly influenced by the snowmelt.This is caused by the high water flux, which results in a high solute flux.To study the effect of the snowmelt on the leaching surface, we also performed a simulation without snowmelt, but with only the precipitation from the year 1997.The weather series is given in Fig. 2b.The high infiltration rate during snowmelt, and the low infiltration rate during winter were removed from the weather series.The day of the solute application is different in Fig. 2a and b, therefore, the infiltration rate is shifted in time.Figure 10a   series 1997 without snowmelt.Leaching starts after 10 mm of drainage, in the highest leaching cells.This is later than in Fig. 9a, because the soil moisture content is higher without the low infiltration rates during winter time with snowmelt.With decreasing leached solute mass (i.e., increasing cumulative area), the cumulative drainage at which solute leaching starts in a cell increases, as was also the case in Figs. 3 and 4. Without snowmelt, the highest solute flux density in a cell is higher than with snowmelt.The solute flux density highly depends on the precipitation rate, leaching only occurs after a precipitation event.In between precipitation events, water may drain, but the solute flux density is much lower than after a precipitation event, as the amount of drainage per day is much lower.
Figure 10b shows the leaching surface of a degradable solute with a pulse, and with the weather series of Fig. 2b.The magnitude of the maximum solute flux density is about 80 times lower as for the tracer (Fig. 10a).This ratio is lower than in the simulations with snowmelt, where the max-imum solute flux density of the degradable solute is similar as of the tracer.The lower solute flux density is caused by the longer residence time in Fig. 10 than in Fig. 9, which leads to more degradation.In contrast to Fig. 9, in Fig. 10, the infiltration rate is low compared to the degradation rate of the solute.Tailing is less for the degradable solute than for the tracer.The tracer still leaches after a precipitation event around 100 mm of drainage, while the degradable solute hardly leaches anymore at that time.The fraction of the soil that contributes to solute leaching is similar for the degradable solute and the tracer, which is about 95 %.

Steady state simulations
As solute leaching is shown to depend on the distribution of the infiltration rate, also a steady state simulation is done with the same random field for the saturated hydraulic conductivity (Fig. 8). Figure 11a shows the leaching surface for a steady state simulation for a tracer.The infiltration rate was 2.5 mm d −1 , which is the average infiltration rate in the snowmelt simulations (averaged over 365 days).In the steady state simulation, the moment of leaching generally increases with decreasing total leaching in a cell, like in the transient simulations.The highest solute flux density is lower than for the transient simulations (0.007 m −2 d −1 ).Approximately 90 % of the cells contribute to the solute leaching, which is less than in the transient simulation during snowmelt (100 %).This is probably caused by the higher soil moisture content during snowmelt in the transient simulations, due to the high infiltration rate.As a result, a larger fraction of the soil is highly conductive.Thus, the leaching is different in a transient simulation than in a steady state simulation.Meyer-Windel et al. (1999) experimentally found that solute breakthrough was similar for transient and steady state conditions in a sandy soil.Kuntz and Grathwohl (2009) found that steady state flow can be used instead of transient flow, except when extreme infiltration events occur.Then, solute leaching was higher in transient simulations than in steady state.On the contrary, in a numerical study Russo et al. (1998) found that transient flow enhances lateral dispersion, mostly at shallow depths.We found that a larger area contributes to solute leaching in the transient simulation, but this is attributed to the higher soil moisture content, not to lateral dispersion.
Figure 11b shows the leaching surface for a steady state simulation for a degradable solute, with an infiltration rate of 2.5 mm d −1 .For the degradable solute, the solute flux density is lower, and the tailing is shorter, due to degradation.The fraction of the soil that contributes to solute leaching is similar for the tracer (88 %) and the degradable solute (86 %) in the steady state simulation.When the marginals of the leaching surfaces in the x-axis are compared, the leaching of the degradable solute is more heterogeneous than of the tracer.This suggests that the heterogeneous soil influences the leaching of the degradable solute more, due to the differences in the travel time, which result in different leached fractions.
To study the influence of the flow rate on the leaching surface, also a steady state simulation with an infiltration rate of 25 mm d −1 was done (Fig. 12a).The solute flux density is higher with an infiltration rate of 25 mm d −1 than of 2.5 mm d −1 , because the water flux is higher.The fraction of the soil that contributes to solute leaching increases with an increasing infiltration rate (96 %). Figure 12b shows the leaching surface of a degradable solute with an infiltration rate of 25 mm d −1 .With a high infiltration rate, the tailing of the degradable solute is more similar to the tailing of the tracer than with a low infiltration rate.When the marginals of the leaching surfaces of the tracer and the degradable solutes in the x-axis are compared, the leaching of the degradable solute is similar as of the tracer, with an infiltration rate of 25 mm d −1 , while they differed when the infiltration rate was 2.5 mm d −1 .With an infiltration rate of 25 mm d −1 , the infiltration rate is high compared to the degradation rate, therefore, there is little time for degradation and thus the spatial distribution of the degradable solute and the tracer are more similar.
In an artificial medium consisting of three different types of sand, Rossi et al. (2008) found that solute mixing between the different types of sand increased with an increasing flow rate.Our results confirm this conclusion.Opposed to a spatially correlated approach, when independent stream tubes are used, the effect of the infiltration rate on the spatial distribution of the solute leaching is ignored, as solutes cannot move laterally.

Effect of depth
In the experiments, the MCS was located at 0.51 m depth.At the field site, the groundwater table is located at 4 m depth.The random field for the scaling factor (Fig. 8) is based on the measurements until 0.51 m.We want to investigate whether, and how, the solute flux density would change with increasing depth.Figure 13 shows the marginal distribution of the x-axis of the leaching surfaces at 0.5 and 1 m depth, both for a tracer and a degradable solute, for the year 1997 with and without snowmelt.Leaching at 1 m depth is more homogeneous than at 0.5 m depth, for cases with and without snowmelt alike, and both for the degradable solute and the tracer.This is in agreement with Persson and Berndtsson (1999), who found that the effect of soil heterogeneity on solute leaching is larger at shallow depth.For snowmelt, the difference between the curves at 0.5 and 1 m depth is smaller than without snowmelt.Thus, with a high infiltration rate, and a high soil moisture content, the depth has little influence on the spatial distribution.With snowmelt, the leaching is more homogeneous than without snowmelt.With the higher infiltration rate, a larger part of the soil is highly conductive.As a result, the leaching is more homogeneous.The same re-sult was found by Persson et al. (2005): in a homogenised soil column the homogeneity in flow increased with an increasing water flux.They stated that a critical soil moisture content might exist.Below this critical content, independent stream tubes might develop, and above it, solute mixing might increase leaching to a more convective-dispersive transport regime.
The spatial distribution for the leaching of the tracer is similar as for the degradable solute, in all cases, except at 0.5 m depth for the simulation without snowmelt.Thus, at smaller depth, the soil heterogeneity is more important for the leaching of degradable solutes.The solute flux density is lower at 1 m depth than at 0.5 m.For the tracer this is caused by dilution and dispersion over a longer distance and time.Besides these effects, degradation lowers the solute flux density for the degradable solute at a larger time.
In the simulations, the properties of the random field for the saturated hydraulic conductivity and the retention curve were uniform.This field was based on the measurements, which were done at 0.51 m depth.Below this depth, the soil heterogeneity might be different than above this depth, attributable to root growth, bioactivity or geology (Pierret et al., 2007;Oades, 1993;French et al., 1994), amongst others.As we have observations until 0.51 m depth, we will assume that the same random field can be used below this depth.

Mean breakthrough and leached mass
From Figs. 9-12 was concluded that the fraction of the soil that contributes to solute leaching increased with an increasing infiltration rate.Here, we will investigate how this conclusion affects the mean breakthrough of a tracer, and the leached mass of a degradable solute.
Figure 14 shows the mean breakthrough time for the model with snowmelt, and without snowmelt, calculated with Eq. ( 4).The points do not follow a 1 : 1 line.mean breakthrough time with snowmelt is lower than without snowmelt, because the infiltration rate during the snowmelt is higher.With snowmelt, the mean breakthrough time increases slower than without snowmelt.This means that with snowmelt a larger fraction of the soil is high conductive than without snowmelt.In the simulation with snowmelt, the solutes were applied with a pulse of six days, while in the simulation without snowmelt the pulse was only one day.When a pulse of one day would be used in the simulation with snowmelt, the effect shown in Fig. 14 is enlarged.The mean breakthrough time, and its distribution influence the leaching of a degradable solute.Figure 15 shows the fraction of the leached mass of a degradable solute divided by the leached mass of a tracer, for each cell.These fractions are compared for the simulation with and without snowmelt.With snowmelt, the fractions are higher, because of the lower mean breakthrough time, as expected.With snowmelt, the fraction decreases slower than without snowmelt.This is caused by the difference in the distribution of the mean breakthrough time, with and without snowmelt.With a stream tube approach this effect would be neglected.As stream tubes are independent and one dimensional, the effect of the infiltration rate on the spatial distribution of the drainage is ignored.This can lead to underestimation of the leaching of degradable solutes.

Conclusions
Two experiments with a multi-compartment sampler (MCS) were done, to investigate the effect of different infiltration rates on the spatial distribution of solute leaching.From the experiment, a random field for a scaling factor for the retention curve was deduced.This random field was used in The standard deviation and the correlation lengths of the random field for the scaling factor can be based on the observations of the experiment.Comparing the spatial distribution of leaching, the model corresponds well with the experiment.The agreement between the observations and the simulations depends both on the standard deviation and the correlation lengths, and on the particular random field.A discretisation of half the correlation length was sufficient to capture all small scale processes.The correlation length in depth did not influence the spatial distribution of the solute leaching much.
The spatial distribution of solute leaching, and which fraction of the soil contributes to solute leaching, is determined by the flow rate.When a stream tube approach would be used, this effect of the infiltration rate would be ignored, as stream tubes are independent.
The infiltration rate largely influences the leaching of the degradable solute.One reason for this is obviously the residence time, which is determined by the infiltration rate.The infiltration rate also determines the fraction of the soil that contributes to solute leaching.Therefore, the infiltration rate also influences the spatial distribution of the solute leaching.For a degradable solute this means that the leaching will be higher than would be estimated with independent stream tubes.
The distribution of the infiltration in time determines the residence time until a control plane at a particular depth.For the case of snowmelt, a steady state simulation with an average infiltration rate would underestimate the leaching of a degradable solute.Without snowmelt, a steady state simulation could overestimate the leaching of a degradable solute.The ratio of the degradation rate over the infiltration rate determines the amount of leaching.

Fig. 1 .
Fig. 1.Experimental setup.Installation of the multi-compartment sampler (MCS), shown from the trench (a).The location of the MCS shown from above (b).The marked square indicates the location of the MCS, 51 cm below soil surface.The roof of the trench is visible in the background.The width of the MCS is 31.5 cm.

Fig. 2 .Fig. 3 .
Fig. 2. Precipitation or infiltration rate (a).1997 with snowmelt, solute application with a pulse of 6 days, snowmelt from day 1 to 6, length of weather series is 367 days (b).1997 without snowmelt, solute application with a pulse of 1 day, length of weather series is 244 days.The lines indicate equal parts of the precipitations series in (a) and (b).(c).snowmelt experiment 2010, solute application with a pulse of 1 day, length of weather series is 62 days.The day is the number of days since the first solute application. 20

Fig. 6 .Fig. 7 .Fig. 6 .
Fig. 6.Sorted observed drainage, and simulated with two standard deviations and correlation lengths for generated random fields with Miller-similarity.

Fig. 15 .Fig. 14 .
Fig. 15.Leached mass of degradable solute (half-life time τ =10 d) divided by the leached mass of the tracer per cell, for the year 1997 with and without snowmelt, at 0.5 m depth.

Table 1 .
Mean and standard deviation of the mean absolute error between observed and simulated sorted and cumulated drainage of five realisations, for the irrigation experiment and the snowmelt experiment.