Similar articles

Abstract: pollutant transport in overland flow over surfaces with spatially varying microtopography, roughness, and infiltration was investigated using the diffusion wave equation and transport rate-based equation. The finite volume method in space and an implicit backward difference scheme in time were employed in the numerical solution of the 2D governing equations. The developed model was first tested against an analytical solution and an experimental study involving overland flow and the associated pollutant transport, subsequently a series of numerical tests were carried out. Non-point source pollution was investigated under spatially varying microtopography, roughness, and infiltration. The simulation results showed that microtopography and roughness were the dominant factors causing significant spatial variations in solute concentration. When the spatially varying microtopography was replaced by a smooth surface, the result was an overestimation of the solute rate at the outlet of the upland. On the other hand, when the spatially varying roughness was replaced by the average roughness and spatially varying infiltration rate by the average infiltration rate, the pollutant discharge at the outlet of the upland was not significantly affected. The numerical results further showed that one cannot ignore the spatial variations of slope and roughness when investigating the local pollutant concentration distribution.

1.  Introduction

 Pollutant in overland flow is a major source of contamination in receiving water bodies and in groundwater due to infiltration. These pollutants arise mostly from surface applied agricultural chemicals, pesticides, and faecal bacteria (Monaghan and Smith, 2012). Solute transfer from soil to surface water due to overland flow may also decrease the efficiency of applied nutrients and chemicals.

 Solute transfer from soil to overland flow and subsequent transport in surface flow were extensively studied both experimentally and numerically (Snyder and Woolhiser, 1985; Wallach et al., 1988; 1989; 2001; Payten and Sanders, 1990; Wallach and Van Genuchten, 1990; Wallach and Shabtai, 1993; Govindaraju, 1996; Garcia-Navarro and Playan, 2000; Yan and Kahawita, 2000; Walton et al., 2000; Abbasi et al., 2003; Deng et al., 2005a; 2005b; He et al., 2009; Ouyang et al., 2012). Most of the previous studies carried out were generally 1D physical experiments in laboratories and 1D solute transport mathematical modeling (Govindaraju, 1996; Garcia-Navarro and Playan, 2000; Yan and Kahawita, 2000; Walton et al., 2000; Wallach et al., 2001; Abbasi et al., 2003). 2D models have also been developed and tested against laboratory experiments and field investigations (VanderKwaak, 1999; Yan and Kahawita, 2007).

 A difficulty that often arises in modeling pollutant transport in the rainfall-runoff process is the determination of the initial value of concentration in overland flow. Many models in previous studies often assume this value to be zero or equal to the initial concentration in soil surface for simplicity and convenience of both the numerical and analytical solutions of pollutant transport, which is not physically reasonable in the pollutant transport processes during rainstorms. In fact, a very high initial concentration caused by the first-rain-flush was observed, instead of the commonly assumed zero value. Deng et al. (2005) developed a transport rate-based model for overland flow and solute transport instead of the commonly assumed zero value.

 In the numerical modeling studies, natural surfaces, over which the transport takes place, are often simplified with smooth surfaces having constant average slope in one direction (in 1D case) or average constant slopes in two directions (in 2D case) (Wallach et al., 1988; 1989; Wallach et al., 1988; Govindaraju, 1996; Singh, 2002). Furthermore, the numerical studies most often assumed homogenous physical properties, thus employing a global roughness coefficient and constant infiltration parameters for the entire field surface (Wallach and Van Genuchten, 1990; Deng et al., 2005a). In the case of planar surfaces with homogeneous properties, such assumptions might be plausible. However, one can rarely find such surfaces in nature. The effects of variability of topography, hydraulic properties, and vegetation density on soil erosion have been studied by modeling and measured data (Tayfur and Singh, 2004; Helmers et al., 2005; Helmers and Eisenhauer, 2006). It is our hypothesis that varying microtopography, roughness, and surface-cover properties of a field may play important roles in the transport of pollutants in overland flow. These characteristics may alter the spatial and temporal variation of pollutants in the field, and into the soil column by infiltration. In addition, these varying characteristics may alter the magnitude of pollutant rates from agricultural fields into the adjacent receiving water bodies. Thus, any modeling studies ignoring the effects of non-planar and non-homogeneous physical effects of the field on the pollutant transport may lead to misleading results. To the knowledge of the authors, so far, there has been no study conducted to investigate effects of varying microtopography, surface-cover, and roughness on pollutant transport in overland flows. Thus, the objective of this study is to numerically investigate effects of varying microtopography, roughness, and infiltration on pollutant transport in overland flows.

2.  Mathematical model

2.1.  2D overland flow model

 The depth-averaged 2D governing equation under the diffusion wave approximation for overland flow can be expressed as (Govindaraju, 1996) (1) , (2) , where h is the water flow depth; R is the rainfall rate; I is the infiltration rate; u and v are the overland flow velocities in x and y directions, respectively; S _fx and S _fy are the friction (energy) slopes in x and y directions, respectively; n_x and n_y are the Manning roughness coefficients in x and y directions, respectively, if anisotropy in roughness is considered; and .

 The diffusion wave model has been shown to provide satisfactory results in most free surface flow applications (Akan and Yen, 1981; Tayfur et al., 1993; Tayfur and Singh, 2004). Infiltration rate in this study was modeled by assuming a constant infiltration rate.

 Solution of Eq. (1) requires initial condition and conditions at the upstream and downstream boundaries. The upstream boundary condition was expressed by the zero flow depth and the downstream boundary condition by the zero flow depth-gradient, which have been commonly used in (Tayfur et al., 1993; Tayfur and Singh, 2004). For the initial condition, a very thin water layer of 1×10⁻⁶ mm for flow depth was assumed to avoid the singularity problem in the numerical solution of Eq. (1).

2.2.  2D pollutant transport model

 During rainstorms, a pollutant in the uppermost layer of the soil is often released to the overland flow. Based on few observations, this study assumes that the solute in this surface active layer are uniformly distributed and the exchange rate of the solute between the active layer and overland flow are proportional to the difference between the concentrations in soil and water volume (Deng et al., 2005b). Generally, pollutant transport in overland flow can be described by the advection-dispersion equation (ADE) (Govindaraju, 1996). However, for overland solute transport, the diffusion/dispersion term can be negligible as compared to other terms (Abbott and Refsgaard, 1996; Wallach et al., 2001). Thus, pollutant transport equation in 2D can be expressed as (3) , where c is the pollutant concentration in overland flow; c _s represents the solute concentration in the active surface layer with a thickness of D _m=εh, with h being the overland flow depth and ε being a dimensionless coefficient; R is the exchange rate of solute between the active surface layer and the overlying flow; J _c is the flux across the water surface per unit surface area; and q _co represents other source or sink terms.

 To avoid the uncertainty in actual initial high level concentration conditions, this study employs the definition of the transport rate of solute, which is the product of surface runoff and solute concentration. By replacing c by c=C/Q _f and c _s by c _s=C _s/Q _f, and splitting each differential term of Ch/Q _f into two terms, i.e., Ch and 1/Q _f, and then multiplying both sides of the equation by Q _f, Eq. (3) yields the following rate-based solute transport equation: (4) , where C is the transport rate of solute by overland flow, C _s is the transport rate of solute from the soil to the overland flow, E=εR, and .

 For a steady and net rainfall flow, in Eq. (4) approaches zero. Decline of the solute transfer rate in the active surface layer caused by the transfer of solutes from soil to the overlying runoff can be assumed to follow an exponential law with a time constant μ as C _s=C ₀exp(−μt) (Deng et al., 2005b).

 The initial condition of Eq. (4) is then given as C(x,0)=0 and C _s(x,0)=C ₀(x), in which C ₀ is the initial transport rate of solute from the soil to the overland flow. The upstream boundary condition was taken as the zero solute transport rate and the downstream boundary condition as the zero gradient of solute transport rate. Eq. (4) can be solved as long as the parameters, E, μ, and C ₀ are determined.

3.  Numerical scheme

 The governing Eqs. (1) and (4) can be rewritten in the following general form: (5) , where Φ represents h and Ch for the surface flow equation and pollutant transport equation, respectively; and S _cΦ includes all source terms in the equation for Φ.

 The equations governing surface flow and pollutant transport are discretised using the finite volume method in space and an implicit scheme in time. In the rectangular 2D control volume, the discretised flow and pollutant transport equations are expressed as (6) .

 Picard iterative technique was applied to handle the nonlinear equations. The discretised equations were solved by the strongly implicit procedure (SIP) of Stone (1968). Accurate numerical solution and mass balance are reached when the convergence criteria of both the Picard iteration and SIP solution loops are satisfied at each time step.

 In the discretized surface equation, variables at the interface of control volume are evaluated by an upstream scheme, which ensures the solution monotonicity. The fluxes in Eq. (6) are functions of the water depth and velocities and are determined with the operator ψ. The operator ψ, which is the nonlinear function of the surface roughness and the magnitude of the water potential (friction slope), is evaluated as a cell parameter at the previous iteration level as (7) , where ε_m is a small number (e.g., 10⁻⁸ to 10⁻¹⁰) that, in the event of a zero fluid potential gradient, avoids a division by zero in Eq. (6) by enforcing a minimum value of Ψ (VanderKwaak, 1999). The pollutant transport is then calculated after flow computations.

4.  Model verification

4.1.  Case 1

 The present pollutant transport model is first tested against an analytical solution for transport of dissolved chemicals in overland flow, given by Rivlin and Wallach (1995). In their mathematical model, the chemical flux, R _ch, across the soil surface was determined as R _ch=k _ch(C _soil−C _water), with k _ch being a mass transfer coefficient, and C _soil and C _water being the concentrations of chemical at the soil surface and the runoff water, respectively (Rivlin and Wallach, 1995; Wallach et al., 2001). The analytical solution of water depth and chemical concentration at the outlet of a 100 m long uniform slope was calculated for a rainfall rate of 1.6 cm/h, a duration of 0.25 h, an infiltration rate of 0.3 cm/h, and k _ch=0.9 cm/h. Fig. 1 shows the comparison of calculated hydrograph and concentration, respectively. A good agreement between the numerical results and analytic solutions for both overland flow depth and chemical concentration at the slope outlet was found with the maximum error less than 8%.

4.2.  Case 2

 The developed model was further tested using an experimental study that investigated overland flow and associated pollutant transport over the soil surface during rainstorms (Deng et al., 2005). The experiment was carried out using a rainfall simulator and a soil flume with a dimension of 3.0 m×0.3 m×0.12 m (length×width×height). Prior to running the experiments with a constant rainfall rate, the soil was wet to field capacity and 0.2 kg of granular salt with a mean diameter of 0.4 mm was applied uniformly onto the soil surface. The saturated water content was 39%, and the initial water content was 20%. During the experiments, surface runoff with solute was collected at the end of the flume. The discharge, salt concentration, and solute transport rate were then measured. The experimental conditions can be found in Table 1. The characteristics of soil and rainfall simulator are given in (de Lima et al., 2003) and the details of experiment can be obtained from Deng et al. (2005a).




 To simulate the experiment, the developed model employed a mesh size of 22×22. The parameters used in the simulation, which were calibrated by Deng et al. (2005b), are listed in Table 2. Fig. 2 shows the comparisons of measured and computed hydrographs and pollutographs. Although the measured discharge varied significantly as shown in Fig. 2a, a relative constant value of 0.055 L/s was obtained. The averaged relative error of measured and simulated discharges was about 5.4%. Figs. 2b and 2c compare the observed and simulated solute transport rates and concentrations, respectively. Both figures show that the simulated results using the developed model are in a good agreement with the experimental data.

5.  Numerical investigation

 To investigate the effects of spatially varying microtopographic surface, roughness, and infiltration rate with their corresponding average values on the pollutant transport dynamics, a set of hypothetical data for these parameters variably distributed over soil flume was generated based on the parameters in Case 2. However, the width of the simulation domain in this section is extended from 0.3 m to 3 m to consider the 2D behavior of pollutant transport. The model simulated the overland flow and pollutant transport over 3 m by 3 m soil flume with a constant slope of 10%, a Manning’s roughness coefficient of 0.025, and a constant infiltration rate of 47.6 mm/h as a basis case.

 The hypothetical spatially varying topographic surface was generated over 3 m by 3 m plane according to the random distribution function with a mean of 0.10 and a standard deviation of 0.029, as shown in Fig. 3a. The range of the generated slope data was 0.05–0.15.




 The hypothetical spatially varying surface roughness was generated over 3 m by 3 m plane according to the random distribution function with a mean of 0.025 and a standard deviation of 0.0089, as shown in Fig. 3b. The range of the generated data was 0.01–0.04, corresponding to bare sand (n=0.010–0.016), gravel surface (n=0.012 to 0.030), bare clay-loam soil (n=0.012 to 0.033), soil with grass (n=0.03 to 0.05), which is physically sound.

 The hypothetical spatially varying infiltration rate was generated over 3 m by 3 m plane according to the random distribution function with a mean of 47.6 mm/h and a standard deviation of 21.33 mm/h, as shown in Fig. 3c. The range of the generated data was 10–85 mm/h, corresponding to loam (I=3.4 mm/h), silt loam (I=6.5 mm/h), sandy loam (I=11 mm/h), loamy sand (I=30 mm/h), and fine sand (I=100 mm/h), which is physically sound.

 Note that one may generate many different random fields of microtopographic surface, surface roughness, and infiltration rate with the same mean and standard deviation. However, to have an insight into the effects of employing varying parameters (as opposed to constant average values) on the pollutant concentration profile, one may employ only one realization, as is discussed in this study.

 Fig. 4 shows velocity, streamline, water depth, solute transport rate, and concentration distributions under the spatially varying microtopographic surface, at a time equal to 128 s when runoff reaches steady. As shown in Figs. 4b and 4c, the concentration and pollutant transport rate have non-smooth profiles, clearly reflecting the impact of the spatially varying slopes which influence the flow depth, velocity, and flow path, as shown in Fig. 4a. Fig. 5 gives solute transport rates and pollutographs at the outlet computed for the cases considering constant and variable slope surface. The transport rates use three sets of variable slopes different in space but with the same mean and standard deviation as shown in Fig. 5a. The variable slopes do not significantly influence the results when they have the same mean and standard deviation. The employment of average constant slopes overestimated the pollutant rate at the outlet by about 9%, while there was no significant influence on concentration.







 Fig. 6 shows water depth, concentration, and solute transport rate under the spatially varying surface roughness, at a time equal to 128 s when runoff reaches steady state. Although both concentration and transport rate profiles are impacted by the variable surface roughness and not smooth, the variable roughness has less impact on pollutant transport than varying slopes. This is because the variation of roughness has less influences on the flow depth, velocity and flow path in comparison with the variation of topography. The concentrations in the upstream are generally higher than those in the downstream, and the solute transport rates show the reverse trend. Fig. 7 compares the concentrations and pollutant transport rates at the outlet under constant and variable surface roughness. The transport rates use three sets of variable roughness different in space but with the same mean and standard deviation (Fig. 7a). The results indicate that the pollutographs computed on the basis of the uniform roughness surface well approximate those obtained from the spatially variable roughness surface.







 Fig. 8 shows water depth, concentration, and solute transport rate under the spatially varying infiltration rate, at a time equal to 128 s when runoff reaches steady state. Similar to the variable roughness, the variable infiltration rate has less impact on both concentration and transport rate profiles, as opposed to the variable slope case. Fig. 9 compares the concentrations and pollutant transport rates at the outlet under constant and variable infiltration rates. The transport rates use three sets of variable infiltration rate different in space but with the same mean and standard deviation (Fig. 9a). The results show insignificant difference between the cases with different infiltration rates.

6.  Conclusions

 A 2D pollutant transport rate-based model combined with the diffusive wave surface flow was developed for simulating overland flow and the associated pollutant transport over surfaces with spatially varying microtopography, roughness, and infiltration. The developed model was tested using a set of analytical solution and laboratory data. The satisfactory simulations of the analytical and experimental data verify the accuracy of the model. The effects of varying slope, roughness, and infiltration on pollutant transport during rainstorms were then numerically investigated by the model. Based on the numerical investigation in this paper, the following conclusions can be drawn:

1. The effects of varying slopes on the flow dynamics, pollutant rates and concentrations are very pronounced. Pollutant load could be overestimated by about 10% if an actual spatially varying microtopography is replaced by constant average slopes. Surfaces are rarely smooth in nature, and the constant-slope assumption is far from reality.

2. The effect of varying roughness on pollutant transport rate and/or concentration at the outlet is relatively negligible, and therefore the use of constant average roughness to represent varying roughness is justified. However, as far as local concentration measures are concerned, the varying roughness should not be replaced by average roughness.

3. Spatially varying infiltration could be represented by spatially averaged infiltration without undue loss of accuracy in estimating pollutant load at the outlet of upland.

---

* Project supported by the National Natural Science Foundation of China (No. 51009120), the Research Fund for the Doctoral Program of Higher Education of China (No. 20090101120065), the State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau of China (No. 10501-243)