We investigate the effect of dynamic boundary conditions on solute transport in unsaturated, heterogeneous, bimodal porous media. Solute transport is studied with two-dimensional numerical flow and transport models for scenarios where either (i) solely infiltration or (ii) more realistic dynamic (infiltration–evaporation) boundary conditions are imposed at the soil surface. Travel times of solute are affected by duration and intensity of infiltration and evaporation events even when cycle-averaged inflow rates of the scenarios are identical. Three main transport mechanisms could be identified based on a criterion for the infiltration rate that is related to the hydraulic conductivity curves of the media. If, based on this criterion, infiltration rates are low, the transport paths for upward and downward transport do not differ significantly, and the breakthrough curves of solute are similar to the one obtained under stationary infiltration. If infiltration rates are moderate, travel paths deviate between upward and downward flow, leading to a trapping of solute and strong tailing of the breakthrough curves. If infiltration and evaporation rates are very high, lateral advective–diffusive transport can lead to very efficient and fast downward transport. Thus, solute breakthrough depends strongly on lateral flow paths enforced by the boundary conditions at the soil surface. If heterogeneity of the materials is not strong and the structure is tortuous, dynamic boundary conditions mainly lead to increased macrodispersion. We test simplified upscaled transport models based on stationary flow rates to estimate breakthrough curves and demonstrate how the transport mechanisms are captured in the model parameters.

Understanding transport of solutes through unsaturated soil in the shallow subsurface is vital to assess nutrient cycling and uptake of water and nutrients by plants, to assess transport of contaminants, to prevent groundwater pollution, and to plan remediation projects (van Dam et al., 2004). The shallow subsurface is characterized by considerable structural heterogeneity of geologic formations that affects the transport of solutes (e.g., Dagan et al., 1990; Jury and Sposito, 1985; Roth, 1995; Roth and Hammel, 1996; Russo et al., 1994; Russo and Fiori, 2008; Russo, 2012; Sudicky, 1986). The coupling to the atmosphere subjects the shallow subsurface to spatially and temporally varying (dynamic) boundary conditions in the form of infiltration and evaporation. As a result of these dynamic boundary conditions, flow in the shallow subsurface is subject to temporal changes in saturation and thus to related structural heterogeneity of the soil (e.g., Roth, 1995; Roth and Hammel, 1996; Kasteel et al., 2009). Additionally, evaporation and infiltration result in temporal changes of flow direction accompanied by possibly differing flow paths during upward and downward flow (Russo et al., 1998). This makes the prediction of transport in the unsaturated zone near the soil surface a challenging problem.

While downward transport of solutes during infiltration events is generally well studied on different scales, upward transport of solutes during evaporation in heterogeneous, artificial porous systems has recently gained attention (Bechtold et al., 2011a, Bechtold et al., 2012; Lehmann and Or, 2009; Nachshon et al., 2011). Heterogeneity in those studies was introduced by two or more sands of different grain size distributions, which are arranged as a composite material with a single material interface or as a stochastic material distribution. Results showed lateral transport and redistribution of solute from coarse to fine sand deeper in the soil column and redistribution from fine to coarse sand at the soil surface. Redistribution of solutes at different depths of the column potentially has an effect on solute breakthrough during subsequent infiltration periods, thus emphasizing the importance of lateral flow for the estimation of travel times. However, infiltration and evaporation are almost exclusively studied independently.

A coherent investigation of the combined effects of infiltration and evaporation is only performed in rare cases (Bechtold et al., 2011a; Russo et al., 1998; Russo and Fiori, 2008; Vanderborght et al., 2006). Russo et al. (1998) numerically simulated three-dimensional field-scale flow and reactive (sorptive) transport in a heterogeneous, unsaturated porous domain to study the effects of soil heterogeneity on solute spreading under realistic conditions. A correlated random Gaussian heterogeneity of the soil properties was induced and infiltration–evaporation periods, as well as water uptake by roots, were considered. Sorption and root water uptake had a retarding effect on solute breakthrough. Compared with steady state conditions, Russo et al. (1998) observed reduced skewing of the solute breakthrough and an increased transversal spread of the solute plume under transient conditions. Furthermore, they noticed that solute spreading in the direction of the mean flow velocity could either increase or decrease because of periodic temporal fluctuations of the mean velocity. Regarding solute breakthrough, an earlier breakthrough, higher peak concentrations and longer tailing was observed. In a more recent numerical study, Vanderborght et al. (2006) confirmed findings of retarded breakthrough and larger transversal redistribution for transient simulations. Lateral solute redistribution from regions with high flow rates to regions with lower flow rates were observed in transient simulations.

In practical applications at field scale, the knowledge of the heterogeneous structures in the shallow subsurface is usually limited. Hence, the direct transfer and representation of natural heterogeneity in numerical models is usually not possible. Furthermore, two- or three-dimensional numerical models that resolve all heterogeneities are computationally demanding. Therefore, prediction of field-scale flow and transport is commonly done by one-dimensional flow or transport models that consider spatially integrated variables rather than exactly representing those (e.g., structural heterogeneity) (Scheibe and Yabusaki, 1998; Vanderborght et al., 2006).

The loss of spatial, structural heterogeneity in these upscaled models is commonly considered by a macrodispersion coefficient or by using double- or multicontinuum approaches adding, for example, a multirate mass transfer term (Li et al., 2011) to the homogeneous model. A thorough discussion about accounting for heterogeneous structure without resolving it can be found in recent literature (e.g., Feyen et al., 1998; Gerke, 2006; Simunek et al., 2003; Vereecken et al., 2007).

Temporal variability of boundary conditions and material properties on the other hand is commonly neglected during the upscaling procedure, which is usually performed assuming steady state flow conditions, for example, by using averaged net fluxes instead of resolving detailed weather data with changing (infiltration–evaporation) boundary conditions. Transient conditions, which have an influence on solute travel times (Russo et al., 1998; Foussereau et al., 2001; Vanderborght et al., 2006; Bechtold et al., 2009), are commonly omitted during the upscaling procedure. However, if upscaled models derived from steady state flow conditions are subsequently used to simulate transient flow and transport, this might lead to erroneous results, depending on the influence of the specific dynamics of the boundary condition.

The objective of this study is to identify and to discuss the impact of dynamic boundary conditions in combination with material heterogeneity that lead to changes of effective solute transport compared with steady state flow conditions. For this purpose, we study test cases with simplified structures using numerical simulations of the coupled flow and transport process. From these cases, we identify three transport phenomena that result from dynamic boundary conditions and that have a strong influence on the effective transport behavior. We test how these phenomena would be captured in macroscopic transport parameters that are based on steady state flow without resolving the fluid distribution. We use the classical advection–diffusion equation and a simple double-continuum model that can be solved semianalytically.

The effect of dynamic boundary conditions is studied using numerical simulations of flow and transport within two-dimensional heterogeneous porous media.

Flow and Transport Equations

Transient flow through the unsaturated porous medium is described by the Richards’ equation (Richards, 1931):
where h is the pressure head [L], being negative in unsaturated conditions, t is the time [T], θ the volumetric water content [L3 L−3], u the pore water velocity vector [L T−1], K the hydraulic conductivity [L T−1], and z the upward-directed vertical coordinate [L]. Volumetric water content, θ(h), and hydraulic conductivity, K(θ), in Eq. [1] are characterized by the functions proposed by van Genuchten (1980):
where θr is the residual, θs is the saturated water content [L3 L−3], α [L−1] and n [-] are parameters of the van Genuchten functions, Ks is the saturated hydraulic conductivity [L T−1], Se [-] is the effective saturation, and l is the tortuosity [-]. Solute transport is described by the advection–dispersion equation (ADE):
where C is the solute concentration in pore water [M L−3] and D is the dispersion tensor [L2 T−1] that is described by,
where αT is the transversal and αL the longitudinal dispersivity [L], Dm is the effective molecular diffusion coefficient [L2 T−1], and I the identity matrix.

Simulation Setup

The cell-centered finite volume flow model Muphi (μϕ, Ippisch et al., 2006) and the random walk particle tracking (RWPT) algorithm PARTRACE (Bechtold et al., 2011b) are coupled and used here to model flow and transport. The ADE (Eq. [4]) is solved with a Lagrangian scheme (RWPT) that, opposed to Eulerian schemes, has proven to give more accurate results for transport in the highly transient shallow subsurface where both advection- and dispersion-dominated transport can occur (Bechtold et al., 2011b, 2012). The scheme’s capability to handle high concentration gradients (e.g., as a result of solute accumulation at the soil surface during upward flow periods or at material interfaces) and its definite mass conservation are especially important in this context.

Two-dimensional saturated and unsaturated flow and transport is simulated in a heterogeneous soil column to study the combined effects of material heterogeneity and time-dependent boundary conditions on travel-time distributions. For this purpose, scenarios that comprise an identical net flux through infiltration and evaporation cycles but differ in intensity and duration of the respective events are studied. Emphasis is put on travel-time distributions at the bottom of the investigated systems that are derived from solute breakthrough curves.

Analogous to recent studies (Lehmann and Or, 2009; Bechtold et al., 2011a), heterogeneity is induced by two different sands (coarse and fine) that are vertically aligned in parallel (Fig. 1a). Thereby, a sharp textural contrast between two hydraulically interacting domains of the porous medium is created. Thus, a preferential flow scenario is set up as it could occur in soils where heterogeneous structure is often subject to strong contrasts, for example, between intra- and interaggregate soil matrix or because of local compaction zones.

This structure is used to have a medium that can clearly be assigned to a double continuum case with preferential flow paths and that is simple enough to allow to identify relevant transport mechanisms. The structure neglects, of course, macroscopic tortuosity referring to the material patterns (which we will refer to as tortuosity in the following), multimodal materials, and other structural features.

To consider tortuosity, a second setting with more complex structure is also analyzed. In this second test case, material heterogeneity is induced by a bimodal structured field (Fig. 1b) that consists of the previously used two sands. Materials in the structured field are arranged such that they remain connected from top to bottom of the domain for some pathways. These connected pathways account for ∼60% of the individual materials’ total area. As a consequence, opposed to the connected pathways, disconnected patches of material also exist. As a result, compared with the vertically layered case, a more complex flow field that enhances overall dispersion is expected.

A shallow system where fluxes in the unsaturated zone are affected by the groundwater table is investigated. The height of the simulated domain is kept constant throughout all simulations at lz = 150 cm. The width of the domain is lx = 20 cm for the first test case with vertically layered structure and lx = 80 cm for the second case where heterogeneity is induced by the structured field. The width is always chosen such that the interaction between the domains, hence the influence of the interface on global flow and transport, is significant. In structures with larger width, effects resulting from this interaction become minor when monitoring space integral values (e.g., breakthrough curves). Boundaries on both sides of the domain are impermeable for flow and transport. For simulations under unsaturated conditions, a groundwater table is set at the bottom of the domain by implementing a Dirichlet boundary condition with h(z = 0) = 0 m. For this unsaturated system, the hydraulic conductivity is not constant for each material but, as a function of water content, variable in space and time. The water retention curve for the applied materials under hydrostatic conditions is obtained by the formula proposed by van Genuchten (1980) (Eq. [2]) and illustrated in Fig. 2a. The associated conductivity curve is illustrated in Fig. 2b.

A spatially uniform flux (Neumann boundary) is prescribed at the top boundary. The typical erratically distributed rain events and evaporation periods are represented by simplified periodic cycles to identify transport mechanisms that result from dynamic change of and lateral flow during upward and downward flow periods. All dynamic simulations comprise an infiltration–evaporation cycle of 10 d (tc) including one infiltration and one succeeding evaporation period. Infiltration duration (ti) is varied from 1 to 96 h (4 d). To obtain a common net infiltration rate, qnet, for all scenarios, the infiltration rate (qi) is calculated as qi = [qnettcqi (tcti)]/ti, thus depending on evaporation rate (qe) and infiltration duration (ti). High evaporation rates (qe) as well as short infiltration duration ti relate to high infiltration fluxes. After prior tests of several net infiltration rates that showed similar behavior, the net infiltration rate of 0.2 cm d−1 is kept constant throughout all simulations in the present work. This means that we cover short and long periods with high and low intensities but all resulting to the same, average (positive) net infiltration flux. A complete overview of the dynamic boundary conditions used in this manuscript can be found in Table 1.

Regarding solute transport, all boundaries are reflective except for the bottom boundary where solute (particles) can leave the domain solely by advective transport. The diffusion coefficient Dm,w in water is set to 2.986 cm2 d−1 and the effective molecular diffusion coefficient for the unsaturated porous medium (Dm) is estimated with the model by Millington and Quirk (1961). Longitudinal dispersivity, αL, is set to 0.1 cm and the transversal dispersivity to αT = 1/10 αL. Diffusion and dispersivity coefficients correspond to commonly used values for sand materials and bromide (e.g., see Bechtold et al. [2012] and Kasteel et al. [2002]).

All simulations had the same initial condition at the time of injection. The initial condition is the equilibrium that established under constant infiltration rate of 0.2 cm d−1. During injection, a constant mass of 1 (dimensionless) is applied as a line source of particles at z = 120 cm above the bottom of the domain. Particles are applied 30 cm below the domain surface to circumvent the effect of redistribution during the first infiltration as a result of large transversal flow in the uppermost part of the domain. At the initial condition, horizontally averaged water contents for the materials at the height of particle injection were 0.055 (coarse) and 0.267 (fine) and thus slightly increased compared with those under hydrostatic conditions (Fig. 3). However, we do not show the model warm up period in the results presented in this manuscript.

The spatial discretization of the grid used for the computation of the flow field is regularly spaced with a nodal distance of 0.5 cm in both horizontal and vertical direction. Figure 1 shows a sketch of the simulated domain including heterogeneity, initial, and boundary conditions. Material properties can be found in Table 2.

To highlight the relevance of changing flow paths on transport, we first outline the case of transport under saturated conditions with changing boundary conditions. The influence of dynamic boundary conditions on cumulative solute breakthrough in heterogeneous media is illustrated exemplarily with the help of three representative scenarios highlighted in Table 1. The goal is to identify transport mechanisms leading to deviations from the transport behavior under stationary infiltration conditions. This will, in turn, allow predicting if and when temporal averaging of boundary conditions (e.g., resulting from natural precipitation time series) is feasible to describe solute breakthrough. Scenarios comprising dynamic infiltration and evaporation in cycles are compared with their stationary infiltration counterpart.

Flow and Transport in Saturated Systems

The vertically layered test case is investigated under fully saturated conditions for two scenarios. The domain is kept fully saturated by applying a large positive pressure head of 10,000 cm to the bottom boundary throughout the simulation. Infiltration duration in these scenarios is varied from 4 d (Scenario 1, Table 1) to 1 h (Scenario 3, Table 1). For comparison, a simulation that comprises solely stationary infiltration with the respective net rate and no evaporation periods is considered. Cumulative breakthrough at the bottom of the domain is presented in Fig. 3. Only minor deviations in the cumulative breakthrough between the stationary infiltration scenario and scenarios comprising dynamic boundary conditions exist, and shape and travel times remain almost identical. Note that the steps displayed in the breakthrough for dynamic scenarios of Fig. 3 are due to the 10-d cycles. The solute breakthrough curve has the typical shape of that of a double-continuum transport model with a clear breakthrough from the fast domain and a tailing resulting from the slow transport out of the slow domain.

Since water content in this system is invariant in time, materials do not experience any temporal or spatial change in hydraulic conductivity. Flow directions and flow paths are simply reversed during downward and upward flux periods. Therefore, solute is transported upward and downward on the same pathways, in this case, preferentially in the coarse material, which has the highest hydraulic conductivity at saturation. In general, except for small deviations, solely temporally averaged net fluxes are relevant for solute breakthrough, and the dynamic boundary conditions play only a minor role.

Crossing of Conductivity Curves in Unsaturated, Heterogeneous Systems

Under steady-state flow conditions, water content and conductivity generally decrease with distance from a groundwater table. The rate of this decrease is determined by the flow rates and by the hydraulic parameters of the materials. Close to the groundwater table, both materials remain close to full saturation. In this region, the coarse material conductivity exceeds the conductivity of the fine material. However, as a result of lower capillary forces, water content and conductivity decrease faster with height in the coarse material than in the fine. Therefore, the conductivity values of fine and coarse material can cross at a certain height above which the fine material has a higher unsaturated hydraulic conductivity than the coarse material (e.g., Roth, 1995; Roth and Hammel, 1996). This crossing point has been experimentally observed at the pedon scale during a tracer experiment under evaporation conditions (Bechtold et al., 2012).

The profiles of hydraulic conductivity have to be considered for the dynamically alternating downward and upward flow. During infiltration periods, starting at the domains surface, the wetting front propagates through the porous medium, causing a local increase of hydraulic conductivity. Evaporation, on the other hand, results in a local drying, causing a local decrease of hydraulic conductivity close to the soil surface. Thereby, conductivity and flow paths become not only variable in space but also in time, leading to potentially changing flow paths during upward or downward flow.

Under realistic conditions, during upward flow, the fine material is the more conductive material close to the soil surface. The coarser material dries out stronger and is thus mostly less conductive than the fine material. The crucial question, whether flow paths close to the soil surface change when the boundary condition switches, thus depends mainly on the preferential flow path during infiltration. If preferential flow during infiltration occurs in the fine material, the flow path remains the same during evaporation. If, on the other hand, preferential flow during infiltration occurs in the coarse material, the flow paths close to the soil surface will change during upward and downward flow.

If a change of flow paths occurs, it can be estimated by comparing the infiltration rate with the hydraulic conductivity curves of the materials. The curves are shown in Fig. 2b and it can be seen that they intersect. At a certain matric potential, the hydraulic conductivities of both materials are equal. For smaller matric potential (wetter soils), the coarse material is more conductive, while for larger (negative) matric potential (dryer soils), the fine material is more conductive. As we hypothesize that the value of hydraulic conductivity at which the two curves intersect is important to predict transport behavior, we use the term intersection conductivity, Kintersect, for it in the following. For steady-state infiltration in a single material and a deep groundwater table, the water content at the soil surface is constant with depth, such that the unsaturated hydraulic conductivity matches the flow rate. If two adjacent layers are infiltrated, the flow rate in the layers adjusts, such that the matric potential in both layers is the same (and constant) and the average of the hydraulic conductivities of both materials is equal to the infiltration flux.

If for the layered case, the infiltration rate would initially be the same in both layers and is higher than the intersection conductivity Kintersect, the matric potential in the coarse material would be lower than in the fine material, such that the fluxes redistribute from the fine to the coarse material, and the preferential flow would be in the coarse material. In this case there would be a change of flow paths during the change from evaporation to infiltration and back.

If, on the other hand, the infiltration rate would initially be constant in both layers and smaller than the intersection conductivity Kintersect, the matric potential in the fine material would be lower than in the coarse material and the fluxes would redistribute from coarse to fine with the preferential path downward occurring in the fine material. In this case, there would be no change of preferential flow paths during upward and downward flow.

To summarize, the preferential flow close to the water table always occurs in the coarse material. Under evaporation conditions, upward preferential flow occurs in the fine material above the depth of conductivity intersection. Under infiltration conditions, however, if the preferential path close to the soil surface changes from fine to coarse or does not depend on whether the infiltration rate is higher or lower than the intersection conductivity.

Flow and Transport in Unsaturated Systems with Double-Continuum Structure

We illustrate the effect of infiltration flow rates quantified by comparing them with the intersection conductivity, Kintersect (which here is 5.8 cm d−1) on transport with the help of three examples chosen from Table 1. The scenarios are chosen such that rates during infiltration periods are below (Scenario 1, mild infiltration), in the vicinity of (Scenario 2, medium infiltration), or above (Scenario 3, strong infiltration) the intersection conductivity Kintersect to illustrate different transport conditions. Consequently, Scenario 1 comprises a scenario with slow infiltration rate for which flow paths during infiltration and evaporation do not change, Scenario 2 is a scenario with medium infiltration rate with changing flow paths; and Scenario 3 is a scenario with fast, strong infiltration rate with strongly changing flow paths and a net upward flow in the fine material (Table 1).

Cumulative breakthrough for Scenarios 1 through 3 under unsaturated conditions is presented in Fig. 4. In general, an earlier first breakthrough for the unsaturated than the saturated system is found for all scenarios and attributed to higher pore (transport) velocity under unsaturated conditions where solute is only transported within the water filled pore volumes of the porous medium. However, it becomes clear that in contrast to the saturated cases, cumulative breakthrough behavior differs significantly between the different scenarios (Fig. 4). Breakthrough of Scenario 1 with mild infiltration is similar to the stationary infiltration case as was found in the saturated example. However, Scenario 2 shows significantly pronounced tailing and Scenario 3 experiences two breakthroughs with less tailing than the second scenario. Deviations between unsaturated scenarios result from temporally and spatially variable water content and associated hydraulic conductivity and change of flow paths. Hence, to understand transport phenomena and deviations of breakthrough between unsaturated scenarios, it is useful to investigate hydraulic conductivity and flow paths, which are illustrated in Fig. 5, at the end of a respective infiltration and evaporation period for each scenario.

Flow Paths in the Different Scenarios

The net infiltration rate of the scenarios is lower than the intersection conductivity Kintersect. In the case of stationary infiltration, the conductivity of the fine material remains thus higher than that of the coarse material close to the soil surface. Therefore, preferential flow in the fine material occurs in this region, which results in a strong lateral flux from coarse to fine material at the soil surface, as a result of spatially constant influx assigned to the top boundary. Preferential flow switches from fine to coarse material in the central region of the domain, which correlates to the intersection point of the hydraulic conductivity curves.

For Scenario 1, comprising a long, mild infiltration with a rate below the intersection conductivity Kintersect, the conductivity of the fine material remains higher than that of the coarse material close to the soil surface during infiltration (Fig. 5a). The flow patterns are thus the same as those for the case of stationary infiltration. During evaporation, flow directions as well as flow paths are reversed (Fig. 5c) analogous to the system under saturated conditions, and in the same way, the breakthrough curve does not differ much from that for stationary infiltration.

With a shorter infiltration duration at a medium infiltration rate in the vicinity of the intersection conductivity Kintersect (Scenario 2), the coarse material, at the end of an infiltration period, throughout the whole vertical extent of the domain, is more conductive (Fig. 5d) than the fine material. Hence, flow redistributes laterally at the surface of the domain and primarily takes part in the coarse material, which prescribes a preferential path, while the fine material is bypassed (Fig. 5e). During evaporation, the preferential path close to the soil surface is in the fine material, and flow paths under upward and downward flow differ from each other.

For a short, very strong infiltration scenario (Scenario 3), where the infiltration rate exceeds the saturated hydraulic conductivity of the fine material (Fig. 5g), the same change of flow paths between infiltration and evaporation as in Scenario 2 occurs. In addition, the wetting front has not propagated through the whole vertical extent of the domain at the end of the infiltration period (Fig. 5h). Therefore, high lateral and downward flow velocities occur mainly above the wetting front close to the soil surface. During evaporation, the component of vertical flow in the fine material is significantly stronger than during infiltration. This leads (in contrast to Scenario 2) to an overall net upward flux in the fine material. As in the second scenario, not only flow direction but also flow paths during evaporation significantly differ from those during infiltration.

To quantify each subdomain’s net contribution to upward and downward flow, vertical pore velocities for the fine material and the coarse material are temporally averaged over an infiltration–evaporation cycle and spatially averaged through the extent of the individual material. The stationary infiltration reference case as well as Scenarios 1 and 2, which comprise weak or medium infiltration events, show a downward directed averaged flow in the fine material with velocities −0.56 cm d−1 (infiltration), −0.45 cm d−1 (Scenario 1), and −0.04 cm d−1 (Scenario 2). Opposed to that, Scenario 3 displays a net upward flow in the fine material with a velocity of 0.16 cm d−1. Cycle- and material-averaged flow in the coarse material is directed downward for all scenarios (infiltration vcoarse = −1.19 cm d−1, Scenario 1 vcoarse = −1.39 cm d−1, Scenario 2 vcoarse = −1.71 cm d−1, Scenario 3 vcoarse = −1.62 cm d−1). Clearly, with increasing strength of infiltration, as long as a steady state regarding downward flux is reached, the coarse subdomains contribution to overall downward flux increases while the fine subdomain contributes more to upward flow.

Solute Distribution and Breakthrough Curves of the Different Scenarios

Solute distribution and main transport directions for Scenario 1 at the end of an evaporation 40 d after the injection of particles is presented in Fig. 6a and after the succeeding infiltration 42 d after injection of particles in Fig. 6b. Solute has moved deeper into the coarse than in the fine soil and shows a higher concentration as a result of strong lateral advective transport from fine to coarse material during infiltration. During evaporation, solute is transported with the mean flow direction and hence upward in the same pathways. This is reflected in the breakthrough curve, which is very similar to that of the stationary infiltration case. The particles are transported mainly along one flow path. This path coincides with the path during stationary infiltration.

For Scenario 2, solute in the coarse material is transported downward fast during infiltration because of preferential flow, while the solute in the fine material is barely affected (Fig. 6c). During evaporation, solute is transported laterally into the fine material and is trapped there (Fig. 6d). For this reason, more solute remains in the fine material for longer time than in Scenario 1. Also, during infiltration, there is no lateral transport from fine to coarse material and hence no fast transport path for particles in the fine material out of the domain. This is reflected in the breakthrough curve. A fast breakthrough of particles occurs earlier than in the stationary infiltration case because of the straight downward path in the coarse material. After the early breakthrough, there is a very strong tailing of the breakthrough curve as a result of the slow release of particles from the fine material. The continuous long tailing, instead of a second distinct breakthrough at later time, indicates that the transport from the fine occurs via exchange with the coarse material.

In Scenario 3 (Fig. 6e,f), solute is transported upward in the fine material over a long distance during evaporation periods and accumulates at the soil surface. Because of a combination of lateral advective transport and diffusion, solute redistributes strongly from fine to coarse material at the soil surface from where it is transported downward fast with the strong infiltration in the coarse material during the next infiltration period. In the breakthrough curve, this results in two distinct breakthrough pulses: the first one coming directly from the particles in the coarse and the second one from particles in the fine material that were transported into the coarse material at the soil surface and washed out one cycle after the first breakthrough pulse.

In summary, it can be stated that upward and downward flow under unsaturated conditions can happen on different flow paths because of the dynamically altering hydraulic conductivity of the subdomains. If and how flow paths change depends on intensity and dynamics of the imposed boundary conditions as well as on the hydraulic conductivity curves of the materials. Pronounced tailing in Scenario 2 can be attributed to trapping of solute in the fine material during evaporation, while the two breakthroughs of Scenario 3 occur as a result of upward transport during evaporation in the fine material followed by lateral transport to the coarse material at the soil surface and fast flushing of the coarse material during infiltration. In conclusion, results from our simulations in bimodal conductivity fields indicate that the use of temporally averaged boundary fluxes is not suitable to describe solute transport in media with a clear double-continuum structure that experience medium or strong infiltration rates (quantified as larger or in the range of the intersection conductivity), such that a change of flow paths during upward and downward flow can be expected.

Predictability of Transport Behavior

An estimation of whether or not breakthrough will be similar to that resulting from stationary infiltration with the same net infiltration rate can be made by the infiltration rate (under dynamic conditions) and the value at which the hydraulic conductivity curves of the materials cross (Kintersect). For infiltration rates lower than that value (our Scenario 1), the behavior will be similar. If the rates are higher (our Scenarios 2 and 3), the behavior will be different. For the higher rates, it is possible that the breakthrough curves exhibit a strong tailing as solute gets trapped in the fine material during evaporation (our Scenario 2). However, if evaporation rates are very high, such that there is a clear net upward flow in the fine material, there could be a very fast breakthrough in few pulses because of the lateral transport at the soil surface (our Scenario 3).

In contrast to Scenario 2, there is no straightforward criterion to predict this behavior. The redistribution process of Scenario 3 requires a net upward flux in the fine material, which has to be large enough to transport particles to the surface. Whether or not this happens is hardly predictable without a numerical flow simulation.

Extension to Additional Scenarios

While the chosen representative scenarios are well suited to illustrate and clearly distinguish between possible occurring transport mechanisms, the evaluation of a wider range of cases (Table 1; Fig. 7) shows that transition between scenarios is continuous. The occurrence of specific transport mechanisms is, besides intensity and dynamic of the boundary conditions, also affected by the spatial distribution and local concentration of particles. Nevertheless, the general behavior follows the mechanisms described above very well. In Fig. 7, the breakthrough curves are plotted for four different evaporation rates. For all cases with infiltration rates clearly below the intersection value (Kintersect = 5.8 cm d−1) of hydraulic conductivities (indicated by bold italic values in Table 1), the breakthrough curve is similar to the one for stationary infiltration. If the infiltration rate is clearly above that value, a preferential downward path occurs in the coarse material. If evaporation is weak (because infiltration time is short), there is a strong tailing of the breakthrough curve. For cases with strong infiltration rates and strong evaporation, there is a redistribution of particles from the fine to the coarse material at the soil surface. These cases usually result in a fast breakthrough and can lead to a second pulse after the first one. To associate the cases in Table 1 with the mechanisms, cases with a net upward flow in the fine material are marked with a double dagger (‡). Not all of these cases show the fast breakthrough, and furthermore, some cases with infiltration rates slightly below Kintersect showed tailing. This leads to the conclusion that the criteria outlined above are ranges rather than sharp criteria for different transport behavior; however, the general trend is confirmed. To quantify the properties of the breakthrough curves, in Table 1 we also give the arrival time for 5% tracer mass (q5). It is clear that the flux boundary conditions at the soil surface play a crucial role for these flow and transport patterns. With Dirichlet boundaries at the top, no lateral flow at the soil surface would take place. However, there is experimental evidence for near-surface lateral flows (Bechtold et al., 2011a). Coupled atmospheric–porous-media two-phase flow simulations are necessary for a more detailed understanding of upper boundary fluxes, their dependency on the heterogeneity of the porous medium, and the lateral flows generated by this interaction. This analysis is beyond the scope of this paper.

Unsaturated Flow and Transport in Bimodal Fields with More Complex Structure

The fields considered above are idealized, as they have two very clearly separated continua with a simple structure. To consider effects resulting from tortuosity, a test case with a tortuous field (Fig. 1b) is implemented. Boundary and initial conditions and material properties are set equal to those used in the vertically layered, unsaturated case. The ratio of fine to coarse material within the domain is kept close to 1. However, the total domain width is increased and set to 80 cm to have structure widths comparable with the vertically layered media test case.

As the parameter contrast of the two materials is not very strong, the tortuosity blurs the clear preferential path structure. A clear upward or downward path is absent, but instead, many tortuous paths occur. Nevertheless, flow paths change for some cases between upward and downward flow in this structure as well. However, the transport in this case can no longer be described by double-continuum behavior but rather by classical advection–macrodispersion behavior. This is seen in the cumulative solute breakthrough at the bottom of the domain (Fig. 8). No strong initial breakthrough or pronounced tailing, which are typical for preferential flow, are existent in the depicted breakthrough curves.

Since lateral flow from coarse to fine material, and vice versa, depends on the retention characteristics and hence material parameters, local redistribution processes in the structured field remain comparable with the vertically layered test case. Lateral flow and solute redistribution from coarse to fine, and vice versa, occur at the same heights as in vertically layered case. Local effects, like trapping of solute in the fine material (Scenario 2) and upward transport in the fine followed by redistribution to the coarse material (Scenario 3), also occur for the structured field. However, because of the system’s strong dispersion, those local transport mechanisms are smeared out. Nevertheless, dynamic boundary fluxes under these conditions still affect solute breakthrough by increasing variances of flow velocity with increasing rates and thus further increasing dispersion (Fig. 8). Although the effect of dynamic boundary conditions is not very strong, it becomes clear that stronger dynamics leads to increased dispersion and earlier breakthrough as a result of the stronger change of flow paths.

By increasing the parameter contrast, the double-continuum behavior with clear preferential paths could also be obtained in the tortuous structure (not shown). In this case, the breakthrough curves and influence of the dynamics look similar to the vertically layered structure but with less pronounced effects.

In the following, we outline upscaled transport models to estimate the breakthrough curves that do not account for transient flow, unsaturated conditions, or heterogeneous soil structure. It will be outlined how the transient flow conditions influence the effective transport parameters.

Effective Model for Double Continua

In this paper, we study the transport in a bimodal medium with a double-continuum structure, where both materials are connected. In principle, traditional dual-continuum flow and transport models (e.g., Simunek and van Genuchten, 2008) could describe such a system. However, our aim is to find an even simpler approach to estimate the solute travel times that is based on cycle-averaged flow rate only.

The simple model is required to reproduce early breakthrough and pronounced tailing as observed in the two-dimensional double-continuum model discussed in the Two-Dimensional Flow and Transport Under Dynamic Boundary Conditions section.

For this purpose we use a double-continuum transport model that considers one-dimensional advective transport with a velocity that is constant in space and time in one continuum and additionally diffusive exchange with a second, stagnant continuum. The advantage of the model is that the breakthrough curve can be calculated semianalytically, such that no numerical simulation is required. To account for diffusion along the transport paths, exchange between continua is modeled as a diffusive process. We used diffusive exchange to keep the number of parameters low (in contrast to a multi-rate, mass transfer model with free rates) and to have smoother exchange than a single rate as in a first-order rate process. The solution of the simple model is obtained by solving the one-dimensional ADE for the continuum with the fastest velocity, which is here assumed to be Domain 1. Neglecting vertical diffusion and considering continuity of mass fluxes and concentrations at the interface between the continua, and furthermore assuming porosity between materials does not differ, one can write the transport equation for the fast continuum, where concentration is horizontally averaged (denoted by angular brackets):
where c [-] is the relative concentration, v [L T−1] the advective relative velocity in continuum 1, and z [L] the vertical coordinate.
The mass flux over the interface between Continuum 1 and 2 is captured as a sink–source term. It is equal to the change of mass in Continuum 2 over time, where ϕi [-] is the volume fractions of layer i. The immobile concentration c2 in equation (Eq. [6]) follows the transport equation:
The horizontal direction is denoted by x. The x direction will not appear in the solution of the problem. Considering the interface conditions, the horizontal average of c2 can be calculated with the method of Green’s function. The solution averaged over the horizontal direction is a function of the mobile concentration c1 and the horizontally averaged Green’s function g, giving the following expression:
The averaged Green’s function for the diffusive exchange has been discussed in several publications (Carrera et al., 1998; Haggerty and Gorelick, 1995). It has been shown by (Carrera et al., 1998) that the diffusive averaged Green’s function can be approximated by a series of exponential functions, such that the model becomes equivalent to a multirate double-porosity model, where the rates are determined by the diffusion coefficient and the width of the layer.
The expansion parameters are
for a slab geometry (Carrera et al., 1998). The tD variable is the diffusive time scale. For a slab with width L this would be tD = L2/Dm. Using Eq. [8] and [9], Eq. [6] can be solved in Laplace space. The Laplace transformed Eq. [6] reads
The angular brackets are here omitted, the tilde denotes the Laplace transformed variable and s is the Laplace variable. This can be solved and yields with c1 (t = 0) = 0 in the domain of interest and for a point injection
With the initial concentration cini = 1, the solution in the stagnant domain is in Laplace space

The solution can numerically be backtransformed with an Iseger algorithm (Iseger, 2009).

The simple model is parameterized by the three parameters v, tD, and ϕ1(with ϕ2 = 1 − ϕ1). The advective velocity v determines the time of the early breakthrough. The time tD represents diffusion and hence has an effect on the shape of the tailing in the breakthrough curve. The volume fraction ϕ1 determines if a larger or smaller part of the transport takes place in the respective continuum, hence leading to a more or less pronounced first breakthrough.

We illustrate the influence of the three model parameters as well as model capabilities and limitations by breakthrough curves in Fig. 9 for the saturated and unsaturated reference cases and representative scenarios. For the vertically layered saturated case (Fig. 9a), no advective cross-transport affects the particles that are injected 30 cm below the surface. Therefore, we assume that both layers contribute equally to the transport and consequently choose ϕ1 = 0.5. We assume that the advective velocity is proportional to the hydraulic conductivity of the coarse material and calculate transport velocities from net infiltration rate as
with the arithmetic mean forumla [L T−1] of both materials. For the diffusive time scale we chose tD = L2/Dm, which is commonly used, with L the width of a layer and Dm = 2.986 cm2 d−1 as used in the two-dimensional numerical simulations of this manuscript. With the above-mentioned parameters, derived from the medium properties, breakthrough calculated by the effective model already matches the breakthrough simulated by the saturated two-dimensional model well without further adjustments or a fitting process (Fig. 9a).
Under unsaturated flow conditions, where flow paths are more complex than in the saturated case, a reasonable estimation of parameters becomes more challenging. With spatially and temporally variable water contents and conductivities, transport velocities are not precisely predictable and preferential flow paths can cross materials inducing solute exchange. However, we are able to estimate parameters for the effective model a priori from material conductivities and inflow rates. As illustrated in the Flow and Transport in Unsaturated Systems with Double-Continuum Structure section, the unsaturated infiltration reference case and Scenario 1, both comprising mild infiltration rates (below the intersection conductivity Kintersect) have a tortuous preferential flow path, which is in the fine material close to the soil surface and moves into the coarse material in the middle of the domain. Because of the twist in the flow path, the effective preferential transport path takes more than half of the domain and thus contributes more to overall transport. Hence, one could assume that a larger portion of advective transport occurs in the fast material. Thus, for the effective model a higher ϕ1 than ϕ1 = 0.5 is chosen. Furthermore, transport velocities for unsaturated cases will increase as a result of lower saturation. The advective transport velocity of the upscaled model is consequently calculated by
where θcoarse is 0.258 and obtained by averaging θ(h) of the water retention curve (Fig. 2). The diffusive time scale was not changed. If ϕ1 is (arbitrarily in the range of 0.5 to 1.0) chosen to 0.8 to compute breakthrough in the effective model for the unsaturated infiltration reference case (Fig. 9b) and the dynamic, mild infiltration Scenario 1 (Fig. 9c), we obtain a good match to the breakthrough of the respective two-dimensional heterogeneous model. Arguably, the ϕ1 value was only estimated to be >0.5, but the exact value could not be estimated. For Scenario 2, fast advective downward transport takes place in the coarse material, where, because of larger infiltration rates (qi), transport velocities are increased. While higher transport velocities could be considered by introducing a scaling factor in the estimation of the upscaled models advective velocity, this would lead to new unclear parameters and is omitted for now, and the same velocity as for the previously mentioned unsaturated scenarios is used. As there is no twist in the flow path, we assume that both materials contribute equally to leaching in this scenario hence selecting ϕ1 = 0.5. The result of the effective model with those parameters is given in Fig. 9d. As suggested the initial breakthrough is not matched well because advective velocities are too slow in the upscaled model.

In Scenario 3, evaporation is so strong that the net flux is upward in the fine domain and particles are thus distributed at the soil surface into the coarse domain. Hence, the preferential transport path is extended. This is captured by an increased volume percentage of the fast domain, ϕ1. The infiltration flux during short, strong infiltrations, and hence transport velocity in the coarse material, increases further. However, we keep v from the previous unsaturated scenarios and solely adjust ϕ1 = 0.8 in the effective model (Fig. 9e). Obviously, the upscaled model is not capable of reproducing the two breakthroughs of the two-dimensional heterogeneous simulations, while the tailing is matched well. Analogous to Scenario 2, the goodness of the initial breakthrough could be matched by increasing the advective velocity in the upscaled model. This would however introduce unknown variables. All of the breakthrough curves in Fig. 9 could be fitted to match the two-dimensional breakthrough curves better; however, the purpose here was to use parameters that are estimated in advance and evaluate the quality of the prediction.

Effective Model for Tortuous Media

For the more tortuous media, the breakthrough curves (Fig. 8) look similar to classical advection–macrodispersion breakthrough curves, and for this reason, we use the classical one-dimensional ADE and solve it analytically with the stationary net infiltration flux divided by porosity for the advective transport velocity. Under fully saturated conditions, the ADE solution with Dm = 2.986 cm2 d−1 and v = qnet/forumla resembles the breakthrough of the two-dimensional model well (not shown). Under unsaturated conditions, higher transport velocities and increased dispersion occur. For unsaturated dynamic scenarios, dispersion and flow velocity further increase proportional to the infiltration rate from Scenarios 1 to 3. The analytical solution of the ADE matches the breakthrough of the two-dimensional model well (Fig. 10). Macrodispersion, even though it clearly increases with inflow rates, could not be predicted quantitatively.

Solute transport under transient conditions has been studied for simple bimodal heterogeneous porous structures with a focus on solute leaching (breakthrough curves). One structure, where setup is similar to setups used in recent studies of evaporation (Lehmann and Or, 2009; Bechtold et al., 2011a), comprises a sharp straight vertical material interface and has a clear double continuum structure where preferential paths develop. A second setup comprises a structured field that considers tortuous pathways. Dynamic boundary conditions were applied at the top of the domain to simulate dynamic atmospheric conditions. The erratic rainfall patterns that occur in nature were simplified to repeating cycles of infiltration and evaporation. Duration and intensity of infiltration and evaporation, phases were varied, while cycle duration and cycle net flux was kept constant to study their effect on solute breakthrough.

For the system comprising a sharp vertical layer, dynamic simulations generally showed a tendency for earlier solute breakthrough and more pronounced tailing than stationary counterparts that comprised solely infiltration with the same net inflow rate. Locally occurring lateral advective transport, which connects paths of fast and slow downward transport, is identified as a significant factor to influence solute breakthrough. In the tortuous structured field, dynamic boundary conditions increased macrodispersion along with transport velocity.

Three general transport mechanisms could be distinguished where occurrence mainly depends on boundary flow rate and material conductivities. If infiltration rates are smaller than the intersection conductivity (Kintersect) of the materials, the dynamic boundary conditions do not have much influence on the breakthrough curve, as flow paths during upward and downward flow do not change. However, flow paths between infiltration and evaporation change if the infiltration rate is larger than the intersection conductivity, and in this case, a strong tailing is obtained, as the finer material acts as a trapping site. If, however, the evaporation rate and time becomes large enough such that the net flow in the fine material is strongly directed upward, solute is transported to the soil surface and laterally redistributed to the coarse material into the preferential downward path where it is flushed out of the domain fast.

A simple one-dimensional upscaled model is considered to estimate breakthrough curves that take advective transport with constant velocity in one continuum and diffusive exchange with another (stagnant) continuum into account. It has been demonstrated that for the saturated system the breakthrough could be predicted from material properties and stationary net inflow rate. Additionally, effects shown in the breakthrough curves could be matched with the model parameters; however, these parameters are difficult to predict quantitatively. For the tortuous medium, a simple advection–diffusion model with constant flow velocity was sufficient to capture the breakthrough curves.

As outlined, the structures considered in this work are strongly simplified, which was necessary to identify the different transport regimes. It was shown that for structures that do not develop clear preferential flow paths, the transport behavior is different and dynamic boundary conditions lead mainly to increased macrodispersion. Hence, to what extent the conclusions drawn here hold for situations when high and low connectivity zones are not connected with each other requires further investigation. It would be necessary to have a prediction for preferential flow behavior based on appropriate measures for heterogeneity to be able to predict the influence of dynamic boundary conditions on transport in more complex and more realistic heterogeneous media. Furthermore, the temporal patterns of infiltration and evaporation used here were strongly simplified periodic patterns. It needs to be tested in future work how the measures used here (rates and time periods) can be transferred to the stochastic measures that are used to describe realistic rainfall patterns.

It has also been pointed out that one essential component of the transport behavior under strong infiltration and evaporation, where solute is transported very fast out of the medium, is the lateral advective transport at the soil surface. This lateral transport is imposed by the constant flux boundary condition at the soil surface. It is an open question how realistic such a boundary condition is. There are indications that this condition is realistic (Bechtold et al., 2011a) as long as evaporation is not reduced because of breaking water films. However, a better understanding of the lateral transport processes directly underneath the soil surface would be needed to assess if the fast transport out of the system is realistic. A coupled surface–subsurface flow and transport model would be needed for this purpose, which also includes an energy balance to capture the evaporation correctly.

This work was conducted in the context of the DFG Project FOR 1083 (NE 824/6–2): MUSIS (Multiscale Interfaces in Unsaturated Soil). Therefore, the authors wish to thank the Deutsche Forschungsgesellschaft (DFG) for funding.

This is an open access article distributed under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).