Next Article in Journal
Hydrodynamic Drivers of Dissolved Oxygen Variability within a Tidal Creek in Myrtle Beach, South Carolina
Previous Article in Journal
Identifying Sensitive Model Parameter Combinations for Uncertainties in Land Surface Process Simulations over the Tibetan Plateau
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Simple Modelling Framework for Shallow Subsurface Water Storage and Flow

1
School of Geosciences, University of Aberdeen, Aberdeen AB24 3UF, UK
2
State Key Laboratory of Soil and Sustainable Agriculture, Institute of Soil Science, Chinese Academy of Sciences, Nanjing 210008, China
3
School of Biological Sciences, University of Aberdeen, Aberdeen AB24 3UU, UK
4
Department of Environment and Geography, University of York, York YO10 5NG, UK
*
Author to whom correspondence should be addressed.
Water 2019, 11(8), 1725; https://doi.org/10.3390/w11081725
Submission received: 9 July 2019 / Revised: 8 August 2019 / Accepted: 14 August 2019 / Published: 19 August 2019
(This article belongs to the Section Hydrology)

Abstract

:
Water storage and flow in shallow subsurface drives runoff generation, vegetation water use and nutrient cycling. Modelling these processes under non-steady state conditions is challenging, particularly in regions like the subtropics that experience extreme wet and dry periods. At the catchment-scale, physically-based equations (e.g., Richards equation) are impractical due to their complexity, while conceptual models typically rely on steady state assumptions not found in daily hydrological dynamics. We addressed this by developing a simple modelling framework for shallow subsurface water dynamics based on physical relationships and a proxy parameter for the fluxes induced by non-unit hydraulic gradients. We demonstrate its applicability for six generic soil textures and for an Acrisol in subtropical China. Results showed that our new approach represents top soil daily fluxes and storage better than, and as fast as, standard conceptual approaches. Moreover, it was less complex and up to two orders of magnitude faster than simulating Richards equation, making it easy to include in existing hydrological models.

1. Introduction

The top few centimetres of soils experience large fluxes associated with cycles in wetting and drying and bioturbation from plant roots and fauna. These fluxes, combined with rainfall, impact and compaction can drastically change pore structure over time [1,2]. As aggregate structures collapse and reform, pathways for water transport and pore water storage are affected [3]. The complexity and dynamics of the pore structure of shallow soils [4] are not adequately considered in existing hydrological models, which could lead to errors in predictions. The potential impact of surface pore structure dynamics on hydrological processes is particularly significant in agricultural regions, where tillage, traffic, carbon depletion and periods of bare soil exacerbate temporal variability. Pore structure dynamics include changes to soil bulk density, pore size distribution and the connectivity of pores [1], but temporal changes in these properties and their impact on hydraulic conductivity are often ignored in hydrological modelling [3].
Complex physically-based relationships between soil water content and water fluxes have long been developed, the most widely used being Richards [5] (1931) equation. Modelling packages, such as HYDRUS [6], include these complex relationships, but are limited by computational time and the flexibility to accommodate external processes. Accommodating a dynamic pore structure, including a thin layer of rapidly changing surface soil, to more accurately model surface fluxes would be unfeasible. External processes like root water uptake have limited flexibility in HYDRUS, with parameterization constrained to pre-implemented Feddes and S-shape models, and computations restricted to 1000 days.
Simpler approaches to simulate sub-surface storage and flow processes exist. Under the assumption of unit gradient, Richards equation becomes the simple Darcy–Buckingham law [7,8]. However, this assumption limits applicability to relatively deep soils and/or for long time-scales where the depth-averaged soil water content dictates the vertical outflux from the soil [9,10]. This precludes the use of such models in shallow soils or in regions, such as the tropics, where surface soil structure dynamics strongly affect influx. In models that are more conceptual, water storage and fluxes are often averaged over a given depth, usually the rooting depth of plants (e.g., [11]). They may also use a more generic undefined depth, such as in the Hydrologiska Byråns Vattenbalansavdelning model, HBV-96 [12]. In HBV-96 and other popular conceptual models, (e.g., the Soil & Water Assessment Tool model, SWAT [13]), vertical fluxes are estimated with an empirical semi-linear reservoir type of behaviour. This typically generates percolation from one layer to another as a function of soil water content and a set of empirical calibrated parameters.
Modelling of soil water in agricultural soils needs to accommodate the shallow top layer of soil that is often cultivated and highly temporally variable. This creates additional limitations to the standard simple depth integrated approach, because this highly managed zone of soil will dominate nutrient cycling and movement within a soil profile [14]. When predicting the export of carbon and nitrogen from soil [15], considerable error could result from failing to account for soil structure dynamics, as changes in water availability can have a dramatic impact on C and N turnover [16] and leaching to groundwater [17]. Salvucci and Entekhabi [18] showed that rainfall events can create highly transient conditions with large soil water fluxes from this upper zone. Similarly, periods of drought can trigger evaporative forcing, large enough to generate temporarily non-negligible upward fluxes across the bottom boundary of a relatively thin top soil zone [19]. In fine-textured soils, this is exacerbated by large capillary forces [20], so assumptions of gravity flow are invalidated. Another challenge to modelling the impacts of surface soil structural properties is their dynamic nature over time due to precipitation [21] or agricultural management practices, such as tillage or the choice of crop [22]. Soil water content in the upper zone is, therefore, a poor predictor of deeper, depth-averaged soil water content, yet this is given inadequate consideration in hydrological modelling. One challenge is the computational requirements of incorporating highly temporal and spatial depth-dependent dynamics.
More simple representations of appropriate shallow soil water storage and fluxes included within existing modelling setups could improve predictions. In this study, we propose such a simple modelling framework of soil water storage and fluxes called the Shallow Subsurface Modelling Framework (SSMF). It accounts for the rapidly changing soil water dynamics in a thin, managed near-surface soil zone under agricultural production, and is parameterised using measurable physical properties. The complexity of the processes it aims at representing is computed in a simple way that directly addresses the physical laws, through a simple mass balance approach, and with a set of three new main features. Since the SSMF was developed with the idea to offer a simple and flexible alternative to the complex solving of Richards equation, we evaluated the performance of the SSMF against Richards equation solutions for a range of soil textures. In order to also anchor the results in real conditions, we also compared the results against field data for a study site in subtropical China. The model’s overall performance, as well as the influence of its features were evaluated.
The drive to develop the SSMF came from field observations at the Sunjia Critical Zone Observatory (CZO) in Jiangxi, China (see, e.g., [23,24]). At this site, large changes in surface soil pore structure and moisture were observed over time, due to seasonal shifts from long dry periods to intense monsoon rainfall. The red clay soils (Acrisols, WRB) at this CZO were derived from kaolinitic quaternary red clay with poor exchange capacity [25]. These soils are highly unstable and prone to erosion [26]. They have pore structures that change depending on vegetation or soil management [27], so there is potential to manage unfavourable conditions and decrease overland flow, erosion and contaminant leaching. Red soils are extremely important globally, covering 20% of each of China, India and Brazil. In China, the red soil area provides food to support 40% of the population [25]. Large inputs of water and fertilizers are typically used to maintain productivity, resulting in considerable leaching to the wider environment [25]. Similar challenges to those observed in Sunjia are faced in other sub-tropical regions. The gradual changes observed in the upper soil structure and their impact on short-term, as well as long-term soil water fluxes and storage are challenging to implement in commonly used models, such as HYDRUS, hence triggering the need to develop a more flexible modelling approach. This study is the first step in the development of SSMF was a proof of concept, and it aimed to validate its basic principles. It involves a detailed explanation of the physics involved in the SSMF structure and equations, and a comparison of its results with those from a Richards equation solution.

2. Site Description and Data

The SSMF framework was developed and tested using hydrological parameters from HYDRUS for six generic soil textures. In addition, an unstable, clay loam soil from south-east China was also used. This soil was from the Sunjia CZO site (116°53′58″–116°54′28″ E, 28°13′45″–28°14′12″ N), which drains an area of 51 ha and is located 4 km northwest of the Red Soil Ecological Experimental Research Station in Yingtan, Jiangxi Province, south-east China (Figure 1). The area has a typically warm, humid subtropical monsoon climate with extreme seasonality. The annual precipitation is ~1800 mm, with a strongly marked rainy season between February and June providing 80% of the annual precipitation. The potential evapotranspiration is ~1200 mm and the mean annual temperature is around 17.8 °C. June, July and August are the warmest months, with 70% of the potential evapotranspiration occurring during this period. The elevation at Sunjia ranges between 41 and 55 m above sea level. Its soil is a “red soil” composed mostly of iron oxides and kaolinite clay, with a clay loam soil texture (referred to hereafter as “CZO clay loam”), following the USDA taxonomy [27]. For the purpose of this study, we focused on the main land use in the Sunjia CZO catchment area during the study period; i.e., upland rainfed fields of peanuts (49.5%). Other key land uses included irrigated paddy fields (20.2%) and citrus trees (17.0%). A more detailed description of the Sunjia CZO site can be obtained from Tahir et al. [23].
All simulations used meteorological data recorded at the Red Soil Experimental Station (116°55′ E and 28°12′ N). Soil moisture data were collected continuously between 23 October, 2012 and 31 December, 2013 using a frequency domain reflectometry (FDR) probes (ML-2x, Delta-T Devices Ltd., Cambridge, UK). These were buried at 0.05, 0.2, 0.4 and 0.8 m depths at three different locations under a peanut crop. The measurement frequency was a 30-min interval, with continuous measurements from 23 October, 2012 to 31 December, 2013, available only at 0.4 and 0.8 m. Soil moisture under the citrus trees has been the subject of study in the past [23], but have not been included in this study for simplicity. Soil moisture was also monitored manually with a time domain reflectometry (TDR) probe (Trime, PICO, IMKO, Ettlingen, Germany) every 15 days from July 2013. These data showed that the soil can be divided into two layers with different hydrological properties; a more permeable layer from 0 to 0.3 m, and another deeper one from 0.3m onwards (analysis not shown). Four runoff plots (6 m by 20 m) were constructed at the up and foot slope positions of citrus and peanut fields, which monitored overland flow (Figure 1). The runoff from these plots was recorded by a tipping bucket system [28] with a counter (Onset Computer Corporation, Bourne, USA).
Modelling extended from 23 October, 2012 to 31 December, 2013, corresponding to soil water content monitoring. In 2013 (January to December 2013), the annual precipitation was slightly less than usual (1584 mm) and the potential evapotranspiration accounted for 60% of the precipitation.

3. Development of the Shallow Subsurface Modelling Framework

3.1. Model Requirements

The SSMF was developed to meet the following requirements:
  • Its results should be comparable to the Richards equation solution;
  • It should be computationally simple and parsimonious, with a minimum number of calibrating parameters;
  • It should be computationally flexible, allowing future accommodation of any type of soil or vegetation process, yet able to compute relatively short-term depth-averaged volumetric soil water content θd (m3/m3) as well as soil water flux Ld (e.g., m/d) across a shallow depth z r (m).

3.2. Model Conceptualization and Development

Figure 2 provides a schematic summary of the general study methodology. The model was tested with the climatic conditions and soil properties of the experimental agricultural study site at the Sunjia CZO. In order to evaluate its viability compared to fully physically-based approaches, and given a known set of soil hydraulic parameters, it was directly evaluated against simulations based on the Richards equation solution from HYDRUS [6] (Figure 2). However, simply calibrating the soil parameters used in the SSMF framework against field data would not allow for the new assumptions to be evaluated independently. As such, the soil hydrological parameters from the van Genuchten soil-water retention function θr, θs, n, l, Ks and α were also obtained from the literature values provided in HYDRUS by the ROSETTA pedotransfer model [29], for six generic soil textures (a loam, a loamy sand, a sandy clay loam, a sandy loam, a silt loam and a silty clay loam), in addition to the CZO clay loam soil. This procedure is described in the Supplementary Material.
The model considers a soil profile of depth z (m) within the unsaturated zone of the soil (Figure 3). In this study the depth extends from the land surface zls (z = 0 m) to a relatively shallow depth, z r (m) so that the dominance of the top surface soil to hydrological behaviour can be adequately represented. Averaging over greater depths, as is common practice in hydrological modelling, would not simulate the strong depth dependent changes in water content observed at the Sunjia CZO with soil sensors.
A range of pedological drivers of hydrological significance could influence the depth of this shallow to surface layer, such as the tillage depth [30], the action of plant roots [31] and disaggregation from intense rainfall. The soil water content considered here was the average volumetric soil water content over the depth z r and the vertical soil water flux was the total water flux across z r . We considered a daily time step. A shorter time step than one day could potentially have been implemented, but in hydrological studies the daily step is commonly used and is relevant to nutrient cycling and agricultural practices.
The first main feature of the SSMF is the introduction of a threshold value of depth-averaged soil water content in the soil layer of interest, θe (m3/m3), above which a negative surface water balance (PE < 0) generates an upward flux across the boundary of the soil layer of interest at depth z r . This process is a simple representation of a non-gravity driven flow process, triggered by a negative water balance at the surface, which may be important for relatively thin layers of soil. If the soil is near saturation in z r , both upwards and downwards water flux across z r occurs, balanced between evaporative forcing and precipitation input into the system. The resulting total flux over a given integration time is therefore the sum of evaporative forcing and precipitation, with each having an opposite sign because of the flux direction. θ e is also indirectly related to the soil water content in the soil below z r , with hydraulic conductivity in both layers being positively correlated at any given time. To restrict the complexity of the SSMF to a minimum, the upwards flux was set to the surface evaporative forcing, and does not account for capillary forces that may play a role. The values of θe are expected to be mostly a function of the soil texture.
The second main feature is the introduction of a proxy parameter a [–] for non-unit hydraulic gradient conditions. A differential in hydraulic pressure between the land surface and z r , which results in a non-unit hydraulic gradient, mitigates the otherwise Darcy flux [7]. Solving the water balance in the soil under non-unit hydraulic conditions is not trivial: It is typically omitted in simple conceptual approaches, limiting their use to conditions under which it can be neglected (deep soil and/or steady-state), or represented explicitly, as in the Richards equation. In the SSMF, the non-unit hydraulic gradient conditions are not accounted for explicitly, but rather approximated with a dimensionless parameter a which mitigates the vertical flux across z r   based on the total value of the previous fluxes. A large total value of previous fluxes relative to the influx at the surface signifies a relatively wet soil profile below z r , and will result in a large negative gradient, thus limiting the amount of water that can be transferred downwards through z r . The parameter a is, therefore, expected to be affected by the soil properties and by the conditions at the surface. In the SSMF, we tested it as a site-specific parameter. Since the unit-hydraulic gradient assumption becomes more valid with large values of a , it is also expected to exhibit a decreasing trend with an increase of z r .
The third main feature of the SSMF, which differentiates it most from existing hydrological models, is a varying computational timestep that is dependent on the relative soil saturation. Conceptually, this results from standard and widely used relationships between soil saturation and hydraulic conductivity, which show that the hydraulic conductivity rapidly increases with soil saturation [32,33,34]. Therefore, periods with high soil water contents can be seen as being in a transient-like state at daily or sub-daily time steps. This becomes even more apparent when the depth of the soil layer is relatively small (i.e., less than a meter) and/or when the soil is fine-textured. In those conditions, hydrological models need to be computed at a very small timestep, at the cost of a long computing time. Total computing time may then become critical for semi-distributed or distributed models, especially at larger (catchment) scales, where a model is required to compute soil water content and vertical fluxes for several dozens to hundreds of plot-scale zones (e.g., [35]). In the SSMF, the computational time step is short when the soil is saturated, reaching a value of d t s a t (hr), and longer when the soil is dry, reaching a value of d t d r y (hr). Between those two values the time step follows a linear relationship. For a daily time-scale, the minimum and maximum values for d t s a t and d t d r y vary between 0 and 24 h, respectively, with d t d r y > d t s a t . The values of d t s a t and d t d r y are expected to be a function of the soil texture and z r .
The routine for the SSMF is as follows. For each day, the potential soil water recharge across the upper boundary of the layer (land surface) was determined from a simple surface water balance q d , p o t l s (m/d), computed from the daily total evapotranspiration E T d (m/d) and the total daily precipitation P d (m/d) as shown in Equation (1):
q d , p o t l s = P d E T d
with
E T d = E T 0 , d . ( θ d 1 θ r θ s θ r ) c
where E T 0 , d is daily total potential evapotranspiration (m/d), θ d 1 (m3/m3) is the soil water content of the previous day, as calculated in the following routine (Equation (8)), θ r is the residual water content (m3/m3), θ s is the saturated water content (m3/m3) and c is an empirical calibrated parameter (no units). This relationship conforms to Laio et al. [36], but does not include a plant transpiration component, where the evapotranspiration is related to the soil moisture thresholds that trigger stomatal closure. In addition, the empirical parameter, c, was added to the original linear equation presented by Laio et al. [36] to gain in generality, and to compensate for the fact that the parameter, Ew, (evapotranspiration at the wilting point) from the original publication is omitted here for simplicity. Rather than using this method to obtain E T d , the SSMF used Penman–Monteith to estimate E T 0 , d .
The daily subroutine was then divided into consecutive time steps t by the following equations.
We first determined how much of q t , p o t l s infiltrated into the soil at the time step t ( q t l s (m/d)). The water that exceeds the available pore space, defined by ( θ s θ t 1 )   z r , for each time step t triggered the computation of a subsequent time step within the daily subroutine (Equation (3)).
q t l s = min ( q t , p o t l s , ( θ s θ t 1 ) . z r ) .
This input of water led to a temporary soil water content θ t , 1 (m3/m3), as defined in Equation (4):
θ t , 1 = ( θ t 1 . z r + q t l s ) / z r
θ t , 1 then determined the time step at which the downwards vertical flux across z r was calculated from the relative soil water saturation according to Equation (4) as:
d t t = min ( t t o t i = 0 i = t 1 d t i , ( d t s a t d t d r y ) . θ t , 1 θ r θ s θ r + d t d r y )
for which d t s a t (hr), and d t d r y (hr), are calibrated parameters, which are the time step at saturation and the time step at a relative saturation of 0, respectively. These parameters are, therefore, unique, time-independent variables calibrated for a given simulation. By definition, d t d r y must be relatively high, while d t s a t must be relatively low, between 0 h and d t d r y . The sum of the subroutine dt values should not lead to a total time of the subroutine greater than the total time of the considered time step ttot (equal to 24 h in the case of the daily time step in this study). The vertical downwards flux across z r for this time step q t , d o w n z r (m/d) was then calculated from the constitutive relationship between soil water content and hydraulic conductivity K t (m/d). In our study, we used the Mualem–van Genuchten relationship to evaluate K t , with the Mualem approximation (m = 1 − 1/n) [37]. For the application to the Sunjia study site, we also considered a modified air-entry value of −0.02 m to represent conditions for the CZO clay loam. Such a small negative air entry value was shown to reduce the effect of non-linearity of the hydraulic conductivity function close to saturated conditions, while still preserving the S-shaped curve typical for the van Genuchten model [38]. This resulted in the relationship presented in Equation (6):
K t ( θ t , 1 ) = K s   ( θ t , 1 θ r θ s θ r ) l [ 1 ( 1 ( θ t , 1 θ r θ s θ r ) 1 / m ) m ] 2
where θ s (m3/m3) is the saturated water content, θ r (m3/m3) is the residual water content, l (no units), and m (no units) are empirical parameters, and K s (m/d) is the saturated hydraulic conductivity.
The vertical downwards flux across z r for the time step t, q t , d o w n z r was then defined as
q t , d o w n z r   = d t t / t t o t . K t . exp ( a . k = 0 k = t 1 q k , d o w n z r d k )
with a [–], a calibrated parameter (here between 0 and 100, Table 1), and, k the dummy variable, used to compute the total previous flux through z r . The term exp ( a . k = 0 k = t 1 q k , d o w n z r d k ) mitigates the vertical flux through zr. This exponential function was chosen to approach the theoretical water-retention relationship between the hydraulic conductivity and the matric potential. Following Equation (7), a low integral value k = 0 k = t 1 q k , d o w n z r d k of the previous vertical fluxes, which indicates very dry conditions in the soil below z r , results in an exponentially-related small vertical flux q t , d o w n z r . In this sense, a can be viewed as a proxy parameter for non-unit hydraulic gradient conditions.
This then allowed for the determination of the depth-averaged soil water content across z r and θ t :
θ t = θ t , 1 q t , d o w n z r / z r
Finally, the potential soil water recharge from the land surface for the next time step q t + 1 , p o t l s was updated
q t + 1 , p o t l s = q t , p o t l s q t l s .
The sub-daily routine stopped when the condition 0 t d t i = t t o t was met. The daily depth-averaged soil-water content value θ d and vertical downwards flux across z r , q d , d o w n z r were, therefore, the values at the last time step of the sub-daily routine (t = tend) and the sum of the downwards flux across z r over the sub-daily routine, respectively, calculated as:
θ d = θ t e n d
q d , d o w n z r = i = 0 i = t e n d q i , d o w n z r
The vertical upwards flux across z r q d , u p z r   (Equations (12a) and (12b)) was triggered when the soil water content was large, above θ e , and when P d E d   < 0:
q d , u p z r =   q d , p o t   l s ,   θ d θ e   a n d   P d E d 0
q d , u p z r =   0 ,   θ d < θ e   o r   P d E d > 0
The total daily resulting vertical flux across z r , L d (m/d) was therefore the sum of the upward flux and the downward flux:
L d =   q d , d o w n z r + q d , u p z r
Eventually when the total time ttot was reached, the excess water constituted the overland flow R d (m/d),
R d =   P d E T d L d Δ S d
where Δ S d (m/d) is the change in the soil water storage over z r ; Δ S d = ( θ d θ d 1 ) / z r .
For the first day of the modelling time period, d = 1, the soil water content was equal to the initial value
θ d = 0 = θ 0
At the first step of the sub-daily routine (t = 0), the value of water recharge across the land surface was the daily potential infiltration q t = 0 , p o t l s = q d , p o t l s when P d E d   ≥ 0. When P d E d   < 0, i.e., the surface evaporative forcing was greater than the precipitation, then the sub-daily routine was computed with water recharge from the surface equal to 0; q t = 0 , p o t l s = 0 . The initial value of the soil water content was the soil water content of the previous day (Equation (16)):
θ t = 0 =   θ d 1

3.3. Model Evaluation

For all seven soils, HYDRUS was used to derive daily soil water content θd,HYDRUS and soil water flux dynamics Ld,HYDRUS (m/d) at eight depths of z r : 0.10, 0.15 0.20, 0.25, 0.30, 0.40, 0.50 and 0.70 m, which were chosen so that the upper soil was more finely represented than the lower soil, where the assumption of a unit gradient is more likely to be valid. Since the CZO clay loam properties are vertically heterogeneous (see Supplementary Material), a weighted arithmetic average as proposed by Destouni [39] was calculated for the calibrated parameters. These averaged soil property parameters were then used as inputs to the SSMF, whereby the initial depth-averaged value over z r of the soil water content was θ 0 = θ d = 1 , H Y D R U S . Next, the SSMF was calibrated against the soil water content θd,HYDRUS through a Monte–Carlo simulation (10,000 runs) and its results Ld,mod were compared with Ld,HYDRUS. Figure 2 provides a summary of this procedure and Table 1 provides the input soil parameters θr, θs, n, l, Ks, α and SSMF parameter calibration ranges for dtdry, dtsat and θ e (as a value relative to the saturated water content θ s , a and c.
Based on the HYDRUS Richards equation simulations, the performance of the SSMF was evaluated using the value of the normalised root mean square deviation (NRMSD), defined as the normalised square root of the variance (no units, Equation (17)).
N R M S D =   d = 1 N x d , S S M F x d , H Y D R U S N / x H Y D R U S ¯
In this relationship N is the total number of days and x refer to the variable of interest.
The results are shown for the best simulation; i.e., the one that minimizes the value of NRMSD for the soil water content time series (in Equation (17), x is then θ); and for the behavioural simulations, defined as the simulations that lead to a NRMSD < 5% of the NRMSD of the best simulation. The model was also evaluated based on the Nash–Sutcliffe efficiency [40] and Kling–Gupta efficiency [41] criteria. The results were very similar to the ones obtained based on the NRMSD, and are, therefore, not presented here.
The resulting time series are presented for the CZO clay loam at z r = 0.40 m which corresponds to one of the depths for which there were continuous soil water content data available at the field site (see Supplementary Material). At z r = 0.40 m the impacts of soil management in the top 20 cm and rapid fluxes due to intense rainfall were greater than the other continuous measurements at z r = 0.80 m. Moreover z r = 0.40 m is very relevant to deep water storage that can be captured by crop roots. The results of the other soil textures and the values of z r in the CZO clay loam are presented in summary figures. As described in Section 3.2, the SSMF is based on the three concepts of: (i) a threshold value of depth-averaged soil water content in the soil layer of interest; (ii) the introduction of a proxy parameter a [–] for non-unit hydraulic gradient conditions; and (iii) a varying computational timestep that is dependent on the relative soil saturation. To evaluate the influence of each of these three main concepts of the SSMF, we evaluate the SSMF’s performance, excluding the parameters one at a time. To evaluate the influence of the parameters a and θ e , they were set to 0 (Equations (7), (12a) and (12b), respectively), and to evaluate the influence of the parameters dtdry and dtsat, the SSMF was computed with a constant hourly timestep (dtt = dt = 1 h, Equation (5)).
The potential evapotranspiration E T 0 , d used in all simulations was calculated using the FAO 56 guidelines [42] based on the daily relative humidity, temperature, wind speed and solar radiation data. To isolate the soil water content dynamics from the effects of the vegetation, the Richards equation solver was calibrated without root water uptake (bare soil conditions) for the CZO clay loam soil. This implies that the soil water fluxes and storage changes are a function of the climate and the soil properties only, which leads to slightly biased soil water retention properties. This is discussed more in Section 4.
The actual evapotranspiration E T d calculated by the SSMF cannot be directly compared to the value obtained from HYDRUS. This is because E T d resulting from Equation (2) is a calibrated value that corresponds solely to the outwards flux from the depth of interest, z r , while the value calculated by HYDRUS is the daily integral over an unknown and time-varying depth. Since the SSMF framework is not designed to be a standalone hydrological model, it can accommodate changes in the representation of the processes that are not the focus of the present study. In particular a different representation of ETd could be implemented, if those fluxes were the focus of interest.

4. Results

4.1. Evaluation of SSMF against Data and HYDRUS Simulations

The calibrated HYDRUS soil water content for the CZO clay loam fell mostly within the range of measured values over the time series (Figure 4b). The value of NRMSD between the average of the measured soil water content time series and the calibrated HYDRUS soil water content was 0.09. The low values of observed data ( θ d < 0.2 ) at the beginning of August were biased by one soil water sensor. Similar to the HYDRUS simulations, the SSMF also mostly fell within the range of the measured θ d (Figure 4b), and the NRMSD between the SSMF values and the data was 0.07, slightly lower than for HYDRUS. The NRMSD between the SSMF values and the HYDRUS values was 0.09.
Although both models simulated the observed data similarly well, the SSMF simulations overestimated the soil moisture measurements during the wetting up phase after the dry season (November 2013). HYDRUS slightly underestimated the soil water content for the first main rainfall events, but its simulations still fell within the range of the measured soil water content values.
For overland flow, the SSMF R d values (Figure 4c) were close to those simulated by HYDRUS (NRMSD < 0.01). Even though the overland flow simulated by the HYDRUS and SMMF models reproduced the major event that occurred on 19 June, neither reproduced any of the small values of observed overland flow prior to that (Figure 4c). There were no overland flow data available from 11 September 2013 onwards. For the vertical fluxes, SSMF L d values were overall underestimated compared to HYDRUS (Figure 4d). Nevertheless, the relative dynamics of vertical fluxes were similar between HYDRUS and the SSMF (Figure 4d, NRMSD < 0.01).
The SSMF was two orders of magnitude faster than the HYDRUS simulations (for example, the simulation for the clay loam soil took 36 s in HYDRUS and less than 1 s with the SSMF on a 2.70 GHz computer), supporting the aim of this study to develop a simplified framework for hydrological modelling that can be applied to complex conditions, such as the depth dependent structural variability and climatic extremes observed in subtropical regions. We also ran the same soil setups in HYDRUS 1D, leading to computation times that were on average two to three times faster than in HYDRUS 2D. However, the analysis was made based on HYDRUS 2D values to allow the validation of the SSMF against a more complex but more realistic representation of the soil water storage and fluxes in the different soil textures.
For each of the seven soil textures and eight different depths, Figure 5 shows the SSMF versus the HYDRUS simulations of soil water content. For the CZO clay loam (CL) down to 0.25 m, θ d values were relatively evenly distributed around the 1:1 line, reflecting an absence of pattern (under/over-estimation). Beyond 0.30 m, the clouds of points reflect a slight skew towards an underestimation of the mid-range values (dots below the 1:1 line). This pattern was, however, less marked for the soil textures. In general, the clouds were typically above the 1:1 line at deeper depths, and mostly for the dry values, reflecting a tendency of the SSMF to overestimate θ d relative to the HDYRUS modelled values. This was most marked for the two silty textures, silty loam, SiL and silty clay loam, SiCL. A narrower cloud centred around the 1:1 line suggests that the wetter values, in all but the two silty textures, were relatively well represented by the SSMF, when compared to the HYDRUS simulations.

4.2. SSMF Parameters Influence

Figure 6 illustrates an extreme case of ignoring θ e or a, or setting dt to a relatively fast 1 h interval to determine impacts on SSMF performance for the CZO clay loam at zr = 0.30 m. Whilst the impact during the wet period was minimal, altering any of these parameters had a significant influence during the long dry period (July to November). During other times, θ d remained high (around 50% higher than when θ e was included in the SSMF) while a resulted in relatively faster and stronger drying of the soil between rainfall events. The SSMF set with a constant dt of 1 h performed similarly regarding the representation of θ d dynamics to when dt was variable in time, within the dtdry; dtsat range. The influences of the parameters on R d   (Figure 4b) and L d (Figure 4d) were not as marked as on θ d , and all setups performed relatively well, with their values and time dynamics remaining close to the values modelled by HYDRUS.
Figure 7 provides a general summary of SSMF performance compared to HYDRUS for all soils, and different depths for each of the SSMF parameters, and the SSMF overall. For   θ d (Figure 7, left column), the SSMF generally performed increasingly better with increasing zr (from NRMSD = 0.09 at zr = 0.10 m to NRMSD = 0.06 at zr = 0.7 m, across all soil profiles), except in the CZO clay loam (CL), where the values of NRMSD slightly increased at 0.30 m. More specifically, where   θ e was set to 0, the NRMSD was consistently poorer, leading to worse performance across all soils and depths (NRMSD = 0.12 on average). If the hydraulic gradient was assumed to be uniform, a = 0 or dt fixed to an 1 h interval, a better NRMSD (NRMSD = 0.08 on average) was observed, almost as good as the complete SSMF (NRMSD = 0.07). Generally, the complete SSMF performed best over individual zr regardless of soil, particularly as zr increased.
For   L d (Figure 7, right column), all the NRMSD values were relatively low (NRMSD < 0.07) except for the silty clay loam (NRMSD > 0.2). The changes in performance with depth were also different from those for   θ d . Unlike for   θ d , there was no clear increasing or decreasing trend in the NRMSD values with the increasing values of zr. For   L d and up to 0.3 m, the setup where a = 0 performed better (NRMSD = 0.03 on average) than the SSMF overall (NRMSD = 0.04 on average). For the deeper depths, the SSMF and the setups where dt fixed to an 1 h interval or with a = 0 performed slightly better (NRMSD = 0.06 on average) than the setup with   θ e   = 0 (NRMSD = 0.07 on average).

5. Discussion

5.1. General Performance of the SSMF

The SSMF reproduced values close to those generated by more complex and computationally intense HYDRUS simulations of both soil water content and vertical fluxes across all values of z r , and in all soil textures (Figure 4, Figure 5 and Figure 7).   θ d values simulated by HYDRUS for dry conditions were slightly higher than those simulated by SSMF. The overall NRMSD between the data and the SSMF was 0.07 while it was 0.09 between HYDRUS and the data, and between HYDRUS and the SSMF for the entire time series across all depths and soil profiles. Wetter soil water contents values were well represented by the SSMF (Figure 5).
In general, the SSMF performed increasingly better with increasing zr. This is consistent with changes in the validity of the underlying assumptions of gravity driven flow and unit hydraulic gradient, which hold better with increasing depth and thickness of a studied soil layer, making the simulations with the SSMF less sensitive towards its specific parameters. That highlights the benefits of the SSMF, particularly in shallow soils. For the CZO clay loam, performance of the SSMF slightly decreased at 0.30 m, and increased again for deeper depths. That reflects the change in soil properties in the CZO clay loam at around 0.3 m, driven by its long-term agricultural management (Table SM2 in the Supplementary Material). Therefore, in the presence of such vertical sharp heterogeneity, the SSMF can be expected to perform less well at said vertical interface. This is a direct result of the proxy parameter a being calibrated with a set of soil hydraulic parameters encountered in above zr and that may not represent accurately the properties lower than zr.
The SSMF was mainly optimized on the soil water content, and not on the vertical fluxes time series, since the soil water content is typically more of interest in most hydrological studies and processes. Nevertheless, the evaporation upwards fluxes of water for all soil textures were mostly well represented, with an average NRMSD of 0.05. However, it can be argued that none of the model efficiency criteria used in this study (NRMSD, Kling-Gupta and Nash-Sutcliffe efficiency criterions) accurately evaluates SSMF vertical flux predictions, as detailed observations were only available for soil water content.

5.2. Evaluation of the Influence of the SSMF Parameters

Optimised performance of the SSMF requires variables θ e and a, with processing speed enhanced considerably by using a variable rather than fixed dt. Of these parameters, θ e when fixed at 0 leads to a large overestimation of the soil water content values, demonstrating its importance to the quality of predictions (Figure 6). This is a key attribute of the SSMF setup. In Equation (12a), if evapotranspiration ( E T d ) exceeds the precipitation ( P d ), as would occur during dry periods, setting θ e to 0 results in an upwards flux, q d , u p z r , that equals the total potential flux ( q d , p o t   l s ). The upwards flux recharging zr is therefore overestimated, so θ d remains relatively high. This impact was greater than using a fixed time step or setting a to 0 (Figure 7), highlighting the value of the θ e parameter.
The parameters dtdry, dtsat, and a had a more marginal influence than θ e on the overall performance for θ d predictions. Their influence was noticeable when daily dynamics were examined. For the CZO clay loam at 0.30 m presented in this study (Figure 6) a representative pattern of strong and fast decreases in soil water contents resulted when a = 0 and dt was set to a 1-h interval. As stated previously, equating the exponential term in Equation (7) to 1 (a = 0), neglects the conditions in the soil below zr (Equation (7)), so the vertical flux equals the hydraulic conductivity corresponding to the conditions in the upper zone zr (unit-hydraulic gradient condition). This causes the differential in soil water content between the upper zone and the soil below zr to be large and positive, producing wetter conditions in the upper soil than in the soil below zr when rain infiltrates through the soil surface. As the simulated vertical flux is overestimated, water subsequently drains too rapidly from the upper zone. A constant value of dt had a similar effect on the water content dynamics as this also resulted in an overestimation of the vertical flux through zr (Equation (7)) for wet water contents.

5.3. Overland Flow

We evaluated the SSMF performance against one of the most commonly used approaches for simulating soil water storage and fluxes; that involved solving the Richards Equation in HYDRUS. By developing a simpler approach that performs as well, we provided a new modelling framework that can be flexibly implemented in a wide range of hydrological models. Nevertheless, evaluating a model against other model simulations assumes that the reference model is right and not biased. In the SSMF, overland flow was implemented in order to close the daily water balance, rather than as an explicit and independent process. In HYDRUS, it is both used to close the water balance, as well as represent it independently as Hortonian flow. Since the overland flow implementations were similar for the two modelling approaches, we assumed that this did not affect our ability to assess how well the other key storage and flow processes were represented. In fact, neither HYDRUS nor SSMF showed significant overland flow. However, overland flow in tropical regions is an important process [43], and leads to large erosion rates at the catchment scale [44]. In the Sunjia catchment, during the wet season and due to the high clay content of the soil, overland flow is generated as a result of infiltration excess [45], which cannot be appropriately accounted for when considering solely a daily water balance model. Therefore, the implementation of an independent overland flow generation process could significantly improve the applicability of the SSMF in subtropical regions. Moreover, the simple binary trigger of upwards flux implemented in the SSMF (with the threshold value of soil water content θ e ) provided insufficient resolution of this process so that the SSMF was not capable of simulating upward fluxes of the smallest magnitudes. This would be important for modelling solute transport, and can be addressed when including the SSMF in a hydrological model that does contain a dedicated overland flow module, such as EROSION 3D [46] or the USDA-Water Erosion Prediction model [47]

5.4. The SSMF Concept Application in the Context of Soil Hydraulics Modelling

The SSMF was developed to meet three key requirements (Section 3.1), derived from current limitations of the modelling of soil hydraulics (Introduction). Firstly, its results should be comparable to Richards equation solution. We showed that the SSMF results were indeed in most cases comparable with the solution to the Richards equation provided by HYDRUS. The optimization based solely on θ d lead to relatively low NRMSD, even for L d . The second requirement involved a parsimonious approach. The SSMF computation was simple, and while the three sets of parameters (dtdry and dtsat, a, and θ e ) together lead in almost all the soil profiles and all values of zr with the best performance, θ e remained the most critical. A relatively short but constant value of dt could suffice to obtain overall good results when the absolute daily values are not of primary importance. Finally, we have provided an approach which is computationally flexible. The SSMF Matlab implementation is provided in the Supplementary Material. It is light and flexible, and its set of equations, presented in this study, can easily be included in any numerical model. The running time is also very fast (two orders of magnitude faster than HYDRUS 2D). In addition, while more simple than other approaches like HYDRUS, the SSMF is parameterised using measurable physical properties, so that their dynamic nature, and the role these play on hydrological processes, although not specifically addressed here, can still be considered.
For the SSMF to be implemented in larger scale models, spatial variability would have to be carefully considered, in particular because SSMF was compared here with, and validated against, point scale measurements and a Richards equation solution. Beven [48] and more recently Or et al. [49] pointed out that Mualem–van Genuchten relationships, and therefore the Richards equation, rely on the assumption of capillary equilibrium with soil water content. This limits its validity to 0.05–0.4 m of spatial extent, as suggested by Vogel and Ippisch [50]. Nevertheless, a rapid computational time achieved by the SSMF lends itself well to distributed or semi-distributed models that bring together small spatial units. When thousands or more units are to be simulated, even fast computers could not make the more parsimonious approach of SSMF, and such a reduction in simulation time, redundant. The heterogeneity of field scale soil water processes could also be better implemented, building on earlier research using randomly distributed hydraulic conductivity to represent field scale soil water processes [51,52].
In the presence of steeper slopes, horizontal fluxes may have to be implemented in the SSMF. The addition of the threshold soil moisture value similar to θ e , above which horizontal flux is triggered, could, for example, be tested. A thorough comparison again with Richards equation solution in 2D would have to then be carried out. However, the three main features of the SSMF would not change, only their calibrated value against a 2D model differ from the 1D version presented in this study.
The results suggested that the SSMF and its three main concepts are representing the soil water content and vertical fluxes well even in cases where the soil is heterogeneous and comprises distinct layers with different hydraulics properties. This is where simple hydrological models normally have strong limitations. The SSMF has been tested in this study with one vertically heterogeneous soil, the CZO clay loam that presents a relatively sharp change in soil properties around 0.3 m. Although   θ d and   L d were represented relatively well, the NRMSD was slightly poorer below 0.3m (Figure 7), where the simulated values of   θ d by the SSMF were lower than the ones simulated by HYDRUS (Figure 5). Testing the SSMF for a larger range of vertically heterogeneous profiles is the next step to validate the SSMF principles. Calibration may still be needed in silty soil textures, where the SSMF performance was poorer in dry conditions (Figure 5).
In summary, the SSMF results can be seen as a proof of concept that the simple representation of non-gravity-driven flow and non-unit hydraulic gradient, as presented in this study, can give a good overall representation of soil water dynamics for a thin upper layer of the soil. Its parsimonious, flexible and faster approach, as opposed to solutions in more complex models like HYDRUS, could provide a wide range of new opportunities. These include cases with large spatial heterogeneity and implementation in larger scale models that typically resort to the assumption of a unit-hydraulic gradient and limit their use to steady-state conditions. Therefore, the SSMF has provided a middle-ground framework that eliminates the key issues of the two extreme approaches normally available.

6. Conclusions

The Shallow Subsurface Modelling Framework (SSMF) proposed here offers a combination of three new concepts for simulating near surface water storage and flow dynamics. It performed as well as a more complex and less flexible Richards equation solution computed by HYDRUS at predicting soil water content and vertical fluxes in soil profiles within a relatively thin surface layer.
We tested this new framework for a data-based clay loam red soil, for a range of eight different soil layer depths (0.10 to 0.70 m). The results were in good agreement with simulations from the widely used, but more complex, Richards equation solution in HYDRUS, and fell within the range of the recorded data. We also tested the SSMF for six homogeneous generic soil profiles of various textures. That also showed relatively good performance for the soil water content and adequate performance for the vertical fluxes. Therefore, the three main features of the proposed approach could be introduced in other hydrological models to represent the processes occurring in the upper soil more accurately than current simple approaches and more efficiently than complex ones. The data necessary for applying this framework are widely available, although a calibration using a soil-water-content time series would be advised. Future improvements may include the computation of a more elaborate relationship between the upwards flux magnitude, soil water content and evaporative forcing. Being computationally simple and flexible, yet with results comparable to the Richards equation solution, the approach has the potential to be implemented in complex modelling frameworks to represent hydrological processes; e.g., in environments with large spatial variability. Since its results are comparable at a site under a monsoonal climate, we have shown that the SSMF can also be implemented under extreme precipitation patterns, such as those found in sub-tropical environments.

Supplementary Materials

The following are available online at https://www.mdpi.com/2073-4441/11/8/1725/s1: S1. HYDRUS 2D set up and calibration. We provide in the Supplementary Material the Matlab code of the developed model (ssmf.m), its subdaily routine (ssmf_subdaily.m) and the pedotransfer equations for different soil textures (pedotransfer_func.m). This is original content, developed by Lucile Verrot, made freely available as part as this study, and developed in Matlab 2014a.

Author Contributions

Conceptualization: L.V., J.G. and P.D.H.; methodology: L.V., J.G. and P.D.H., coding: L.V.; validation: all authors; formal analysis: L.V., J.G. and P.D.H.; resources: P.D.H., G.Z. and X.P.; data curation: G.Z., L.G. and X.P.; writing—original draft preparation: L.V. and J.G.; writing—review and editing: all authors; project administration: P.D.H., G.Z., J.G., J.U.S. and X.P.; funding acquisition: P.D.H., G.Z., M.E.H., X.P., J.U.S. and J.G.

Funding

This study was part of the UK–China Red Soils CZO project, funded by the National Environment Research Council (grant NE/N007611/1) and the National Sciences Foundation of China (NSFC: 41571130051, 41571130053, 41371235).

Acknowledgments

Special thanks go to the staff of the Ecological Experimental Station of Red Soil, of the Institute of Soil Sci. of CAS, who provided the detailed meteorological data. We would like to thank three anonymous reviewers for providing constructive criticism on an earlier version of this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Augeard, B.; Bresson, L.M.; Assouline, S.; Kao, C.; Vauclin, M. Dynamics of soil surface bulk density: Role of water table elevation and rainfall duration. Soil Sci. Soc. Am. J. 2008, 72, 412–423. [Google Scholar] [CrossRef]
  2. Liu, H.; Lei, T.W.; Zhao, J.; Yuan, C.P.; Fan, Y.T.; Qu, L.Q. Effects of rainfall intensity and antecedent soil water content on soil infiltrability under rainfall conditions using the run off-on-out method. J. Hydrol. 2011, 396, 24–32. [Google Scholar] [CrossRef]
  3. Ahuja, L.R.; Ma, L.W.; Timlin, D.J. Trans-disciplinary soil physics research critical to synthesis and modeling of agricultural systems. Soil Sci. Soc. Am. J. 2006, 70, 311–326. [Google Scholar] [CrossRef]
  4. Schwen, A.; Zimmermann, M.; Bodner, G. Vertical variations of soil hydraulic properties within two soil profiles and its relevance for soil water simulations. J. Hydrol. 2014, 516, 169–181. [Google Scholar] [CrossRef]
  5. Richards, L.A. Capillary conduction of liquids in porous media. Jpn. J. Appl. Phys. 1931, 1, 318–333. [Google Scholar] [CrossRef]
  6. Šimůnek, J.; van Genuchten, M.T. Using the Hydrus-1D and Hydrus-2D codes for estimating unsaturated soil hydraulic and solute transport parameters. In Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media; van Genuchten, M.T., Leij, F.J., Wu, L., Eds.; University of California: Riverside, CA, USA, 1999; pp. 1523–1536. [Google Scholar]
  7. Buckingham, E. Studies on the Movement of Soil Moisture; Bull. USDA, Bureau of Soils: Washington, DC, USA, 1907; p. 38.
  8. Darcy, H. Détermination des lois d’écoulement de l’eau à travers le sable. In Les Fontaines Publiques de la Ville de Dijon; Victor Dalmont: Paris, France, 1856; pp. 590–594. [Google Scholar]
  9. Gandolfi, C.; Facchi, A.; Maggi, D. Comparison of 1D models of water flow in unsaturated soils. Environ. Model. Softw. 2006, 21, 1759–1764. [Google Scholar] [CrossRef]
  10. Verrot, L.; Destouni, G. Data-model comparison of temporal variability in long-term time series of large-scale soil moisture. J. Geophys. Res.-Atmos. 2016, 121. [Google Scholar] [CrossRef]
  11. Porporato, A.; D’odorico, P.; Laio, F.; Ridolfi, L.; Rodriguez-Iturbe, I. Ecohydrology of water-controlled ecosystems. Adv. Water Resour. 2002, 25, 1335–1348. [Google Scholar] [CrossRef]
  12. Bergström, S. The HBV model. In Computer Models of Catchment Hydrology; Singh, V.P., Ed.; Water Resources Publications: Highlands Ranch, CO, USA, 1995; pp. 443–476. [Google Scholar]
  13. Santhi, C.; Srinivasan, R.; Arnold, J.G.; Williams, J.R. A modelling approach to evaluate the impacts of water quality management plans implemented in a catchment in Texas. Environ. Modell. Softw. 2006, 21, 1141–1157. [Google Scholar] [CrossRef]
  14. Deng, J.; Zhu, B.; Zhou, Z.; Zheng, X.; Li, C.; Wang, T.; Tang, J. Modeling nitrogen loadings from agricultural soils in southwest China with modified DNDC. J. Geophys. Res.-Biogeo. 2011, 116. [Google Scholar] [CrossRef] [Green Version]
  15. Dawson, J.J.C.; Smith, P. Carbon losses from soil and its consequences for land-use management. Sci. Total Environ. 2007, 382, 165–190. [Google Scholar] [CrossRef] [PubMed]
  16. Manzoni, S.; Taylor, P.; Richter, A.; Porporato, A.; Ågren, G.I. Environmental and stoichiometric controls on microbial carbon-use efficiency in soils. New Phytol. 2012, 196, 79–91. [Google Scholar] [CrossRef] [PubMed]
  17. Zhu, Z.L.; Wen, Q.X.; Freney, J. (Eds.) Nitrogen in Soils of China; Springer Science & Business Media: Berlin, Germany, 2012; Volume 74. [Google Scholar]
  18. Salvucci, G.D.; Entekhabi, D. Equivalent steady soil moisture profile and the time compression approximation in water balance modeling. Water. Resour. Res. 1994, 30, 2737–2749. [Google Scholar] [CrossRef]
  19. Brutsaert, W. The daily mean zero-flux plane during soil-controlled evaporation: A Green’s function approach. Water. Resour. Res. 2014, 50, 9405–9413. [Google Scholar] [CrossRef]
  20. Lehman, P.; Or, D. Evaporation and capillary coupling across vertical textural contrasts in porous media. Phys. Rev. E 2009, 80, 46318. [Google Scholar] [CrossRef] [PubMed]
  21. Peng, X.; Horn, R.; Smucker, A. Pore shrinkage dependency of inorganic and organic soils on wetting and drying cycles. Soil. Sci. Soc. Am. J. 2007, 71, 1095–1104. [Google Scholar] [CrossRef]
  22. Or, D.; Ghezzehei, T.A. Modeling post-tillage soil structural dynamics: A review. Soil Till. Res. 2002, 64, 41–59. [Google Scholar] [CrossRef]
  23. Tahir, M.; Lv, Y.; Gao, L.; Hallett, P.D.; Peng, X. Soil water dynamics and availability for citrus and peanut along a hillslope at the Sunjia Red Soil Critical Zone Observatory (CZO). Soil Till. Res. 2016, 163, 110–118. [Google Scholar] [CrossRef] [Green Version]
  24. Lv, Y.; Gao, L.; Geris, J.; Verrot, L.; Peng, X. Assessment of water sources and their contributions to streamflow by endmember mixing analysis in a subtropical mixed agricultural catchment. Agr. Water. Manag. 2018, 203, 411–422. [Google Scholar] [CrossRef]
  25. He, Z.; Zhang, M.; Wilson, M.J. Distribution and classification of red soils in China. In The Red Soils of China; Wilson, M.J., He, Z., Yang, X., Eds.; Springer: Dordrecht, The Netherlands, 2004; pp. 29–33. [Google Scholar]
  26. Ma, L.; Cheng, Y.; Wang, J.; Yan, X. Mechanical insights into the effect of fluctuation in soil moisture on nitrous oxide emissions from paddy soil. Paddy Water Environ. 2017, 2, 359–369. [Google Scholar] [CrossRef]
  27. Zhou, H.; Peng, X.; Peth, S.; Xiao, T.Q. Effects of vegetation restoration on soil aggregate microstructure quantified with synchrotron-based micro-computed tomography. Soil Tillage Res. 2012, 124, 17–23. [Google Scholar] [CrossRef]
  28. Soil Science Division Staff. Soil Survey Manual; Ditzler, C., Scheffe, K., Monger, H.C., Eds.; USDA Handbook 18; Government Printing Office: Washington, DC, USA, 2017.
  29. Schaap, M.G.; Leij, F.J.; van Genuchten, M.T. ROSETTA: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 2001, 251, 163–176. [Google Scholar] [CrossRef]
  30. Salem, H.M.; Valero, C.; Muñoz, M.A.; Rodríguez, M.G.; Silva, L.L. Short-term effects of four tillage practices on soil physical properties, soil water potential, and maize yield. Geoderma 2015, 237, 60–70. [Google Scholar] [CrossRef]
  31. Laio, F.; D’Odorico, P.; Ridolfi, L. An analytical model to relate the vertical root distribution to climate and soil properties. Geophys. Res. Lett. 2006, 33. [Google Scholar] [CrossRef]
  32. Brooks, R.; Corey, T. Hydraulic properties of porous media. Hydrology Papers; Colorado State University: Fort Collins, CO, USA, 1964; 24p. [Google Scholar]
  33. Campbell, G.S. A simple method for determining unsaturated conductivity from moisture retention data. Soil Sci. 1974, 117, 311–314. [Google Scholar] [CrossRef]
  34. Van Genuchten, M.T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil. Sci. Soc. Am. J. 1980, 44, 892–898. [Google Scholar] [CrossRef]
  35. Vázquez, R.F.; Beven, K.; Feyen, J. GLUE based assessment on the overall predictions of a MIKE SHE application. Water Resour. Manag. 2009, 23, 1325–1349. [Google Scholar] [CrossRef]
  36. Laio, F.; Porporato, A.; Ridolfi, L.; Rodriguez-Iturbe, I. Plants in water-controlled ecosystems: Active role in hydrologic processes and response to water stress: II. Probabilistic soil moisture dynamics. Adv. Water Resour. 2001, 24, 707–723. [Google Scholar] [CrossRef]
  37. Mualem, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 1976, 12, 513–522. [Google Scholar] [CrossRef] [Green Version]
  38. Vogel, T.; van Genuchten, M.T.; Cislerova, M. Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions. Adv. Water Resour. 2001, 24, 133–144. [Google Scholar] [CrossRef]
  39. Destouni, G. Applicability of the steady state flow assumption for solute advection in field soils. Water Resour. Res. 1991, 27, 2129–2140. [Google Scholar] [CrossRef]
  40. Nash, J.E.; Sutcliffe, J.V. River flow forecasting through conceptual models part I—A discussion of principles. J. Hydrol. 1970, 10, 282–290. [Google Scholar] [CrossRef]
  41. Gupta, H.V.; Kling, H.; Yilmaz, K.K.; Martinez, G.F. Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling. J. Hydrol. 2009, 377, 80–91. [Google Scholar] [CrossRef] [Green Version]
  42. Allen, R.G.; Pereira, L.S.; Raes, D.; Smith, M. Crop evapotranspiration: Guidelines for computing crop requirements. In Irrigation and Drainage Paper; No. 56; FAO: Rome, Italy, 1998; 300p. [Google Scholar]
  43. Zhang, M.K.; Xu, J.M. Restoration of surface soil fertility of an eroded red soil in southern China. Soil Tillage Res. 2005, 80, 13–21. [Google Scholar] [CrossRef]
  44. Zhang, W.L.; Wu, S.X.; Ji, H.J.; Kolbe, H. Estimation of agricultural non-point source pollution in China and the alleviating strategies I. Estimation of agricultural non-point source pollution in China in early 21 century. Sci. Agric. Sin. 2004, 37, 1008–1017. [Google Scholar]
  45. Wei, L.; Zhang, B.; Wang, M. Effects of antecedent soil moisture on runoff and soil erosion in alley cropping systems. Agr. Water. Manag. 2007, 94, 54–62. [Google Scholar] [CrossRef]
  46. Schmidt, J.; Werner, M.V.; Michael, A. Application of the EROSION 3D model to the CATSOP watershed, The Netherlands. Catena 1999, 37, 449–456. [Google Scholar] [CrossRef]
  47. Flanagan, D.C.; Nearing, M.A. USDA-Water Erosion Prediction Project: Hillslope profile and watershed model documentation. In Nserl Report; USDA-ARS National Soil Erosion Research Laboratory: West Lafayette, IN, USA, 1995. [Google Scholar]
  48. Beven, K. How far can we go in distributed hydrological modelling? Hydrol. Earth. Syst. Sci. Discuss. 2001, 5, 1–12. [Google Scholar] [CrossRef]
  49. Or, D.; Lehmann, P.; Assouline, S. Natural length scales define the range of applicability of the Richards’ equation for capillary flows. Water Resour. Res. 2015, 51, 7130–7144. [Google Scholar] [CrossRef]
  50. Vogel, H.-J.; Ippisch, O. Estimation of a critical spatial discretization limit for solving Richards’ equation at large scales. Vadose Zone J. 2008, 7, 112–114. [Google Scholar] [CrossRef]
  51. Bresler, E.; Dagan, G. Unsaturated flow in spatially variable fields. 2. Application of water flow models to various fields. Water Resour. Res. 1983, 19, 421–428. [Google Scholar] [CrossRef]
  52. Dagan, G.; Bresler, E. Unsaturated flow in spatially variable fields. 1. Derivation of models of infiltration and redistribution. Water Resour. Res. 1983, 19, 413–420. [Google Scholar] [CrossRef]
Figure 1. Map of the Sunjia catchment and its land uses; the altitude lines are in 1 m increments.
Figure 1. Map of the Sunjia catchment and its land uses; the altitude lines are in 1 m increments.
Water 11 01725 g001
Figure 2. Schematic representation of the general study methodology; the goal of the Shallow Subsurface Modelling Framework (SSMF) is to use widely available input data of soil water content θd; i.e., soil water retention properties and climate variables (precipitation Pd and evapotranspiration ETd), to predict future soil water content, θd,mod, and the dynamics of water fluxes, Ld,SSMF. To assess the performance of the SSMF, we compared the model outputs with Richards equation solution values, θd,HYDRUS and Ld,HYDRUS.
Figure 2. Schematic representation of the general study methodology; the goal of the Shallow Subsurface Modelling Framework (SSMF) is to use widely available input data of soil water content θd; i.e., soil water retention properties and climate variables (precipitation Pd and evapotranspiration ETd), to predict future soil water content, θd,mod, and the dynamics of water fluxes, Ld,SSMF. To assess the performance of the SSMF, we compared the model outputs with Richards equation solution values, θd,HYDRUS and Ld,HYDRUS.
Water 11 01725 g002
Figure 3. Model conceptualization of the Shallow Subsurface Modelling Framework (SSMF). The fluxes are P, the precipitation at the land surface, E, the evaporation from the soil, R the overland flow and L, the vertical losses through the bottom boundary of z r , the upper zone of the soil. θ is the depth-averaged volumetric water content over z r and SHP stands for Soil Hydraulic Properties over z r . The vertical axis is positively oriented from the surface down the soil profile. The values in blue are the inputs to the SSMF; the values in green are the outputs.
Figure 3. Model conceptualization of the Shallow Subsurface Modelling Framework (SSMF). The fluxes are P, the precipitation at the land surface, E, the evaporation from the soil, R the overland flow and L, the vertical losses through the bottom boundary of z r , the upper zone of the soil. θ is the depth-averaged volumetric water content over z r and SHP stands for Soil Hydraulic Properties over z r . The vertical axis is positively oriented from the surface down the soil profile. The values in blue are the inputs to the SSMF; the values in green are the outputs.
Water 11 01725 g003
Figure 4. Daily observed precipitation (a), soil water content θ d (b), overland flow R d (c), and vertical flux L d (d) in the CZO clay loam profile of the Sunjia catchment at zr = 0.40 m, from November 2012 until December 2013. In (bd), the time series are obtained from the calibrated HYDRUS software (red) and from the Shallow Subsurface Modelling Framework (SSMF, blue). For θ d (b) and R d (c), the data range and its average value across the 12 sensors are represented respectively by a green shaded area and a green line.
Figure 4. Daily observed precipitation (a), soil water content θ d (b), overland flow R d (c), and vertical flux L d (d) in the CZO clay loam profile of the Sunjia catchment at zr = 0.40 m, from November 2012 until December 2013. In (bd), the time series are obtained from the calibrated HYDRUS software (red) and from the Shallow Subsurface Modelling Framework (SSMF, blue). For θ d (b) and R d (c), the data range and its average value across the 12 sensors are represented respectively by a green shaded area and a green line.
Water 11 01725 g004
Figure 5. Scatter plots of the soil water content θ d values [m3/m3] from the Shallow Subsurface Modelling Framework (SSMF) (y-axis) and from HYDRUS (x-axis), across the eight studied depths z r (rows), and for the seven soil profiles (columns): CZO clay loam (CL), loam (L), loamy sand (LSa), sandy clay loam (SaCL), sandy loam (SaL), silt loam (SiL) and silty clay loam (SiCL). The black line represents the 1:1 line.
Figure 5. Scatter plots of the soil water content θ d values [m3/m3] from the Shallow Subsurface Modelling Framework (SSMF) (y-axis) and from HYDRUS (x-axis), across the eight studied depths z r (rows), and for the seven soil profiles (columns): CZO clay loam (CL), loam (L), loamy sand (LSa), sandy clay loam (SaCL), sandy loam (SaL), silt loam (SiL) and silty clay loam (SiCL). The black line represents the 1:1 line.
Water 11 01725 g005
Figure 6. Influence of ignoring θ e or a, or setting dt to a relatively fast 1 h interval on SSMF performance versus HYDRUS. The outputs are (a) daily observed precipitation, (b) overland flow R d , (c) soil water content θ d and (d) vertical flux L d for the CZO clay at zr = 0.30 m. In (bd), the time series are obtained from the calibrated HYDRUS software (red) and from the Shallow Subsurface Modelling Framework (SSMF, blue), as well as for the SSMF, for which one of those parameters was ignored: θ e (pink line), a (black line), and d t d r y and d t s a t (green line).
Figure 6. Influence of ignoring θ e or a, or setting dt to a relatively fast 1 h interval on SSMF performance versus HYDRUS. The outputs are (a) daily observed precipitation, (b) overland flow R d , (c) soil water content θ d and (d) vertical flux L d for the CZO clay at zr = 0.30 m. In (bd), the time series are obtained from the calibrated HYDRUS software (red) and from the Shallow Subsurface Modelling Framework (SSMF, blue), as well as for the SSMF, for which one of those parameters was ignored: θ e (pink line), a (black line), and d t d r y and d t s a t (green line).
Water 11 01725 g006
Figure 7. Normalised root mean square deviation (NRMSD) between SSMF and HYDRUS for the daily soil water content θ d (left column) and for the daily vertical flux L d (right column). Simulations are for the eight soil texture profiles: CZO clay loam (CL), loam (L), loamy sand (LSa), sandy clay loam (SaCL), sandy loam (SaL), silt loam (SiL) and silty clay loam (SiCL). The values of the studied depths (zr) are on the x-axis; the NRMSD values on the y-axis. The NRMSD values are presented for the simulations from the SSMF (blue crosses), and for the setups where the individual parameters were omitted: θ e = 0 (pink circles), a = 0 (black squares) and dt = 1 h (green triangles).
Figure 7. Normalised root mean square deviation (NRMSD) between SSMF and HYDRUS for the daily soil water content θ d (left column) and for the daily vertical flux L d (right column). Simulations are for the eight soil texture profiles: CZO clay loam (CL), loam (L), loamy sand (LSa), sandy clay loam (SaCL), sandy loam (SaL), silt loam (SiL) and silty clay loam (SiCL). The values of the studied depths (zr) are on the x-axis; the NRMSD values on the y-axis. The NRMSD values are presented for the simulations from the SSMF (blue crosses), and for the setups where the individual parameters were omitted: θ e = 0 (pink circles), a = 0 (black squares) and dt = 1 h (green triangles).
Water 11 01725 g007
Table 1. Input soil parameter values and the range of the calibration parameters used in the Shallow Subsurface Modelling Framework (SSMF). The procedure to determine the value of the soil parameters is described in the Supplementary Material. For the CZO clay loam profile, the soil parameters θr, θs, n, l, Ks and α reported in this table were obtained from the calibration of the Richards equation with HYDRUS on soil data from the Sunjia catchment. For the six generic soil texture profiles, they are the literature values provided in HYDRUS by the ROSETTA pedotransfer model [29]. For all soil profiles, dtdry, dtsat a, θe and c are the SSMF parameters. The ranges that were used for the calibration of the SSMF are given in the table.
Table 1. Input soil parameter values and the range of the calibration parameters used in the Shallow Subsurface Modelling Framework (SSMF). The procedure to determine the value of the soil parameters is described in the Supplementary Material. For the CZO clay loam profile, the soil parameters θr, θs, n, l, Ks and α reported in this table were obtained from the calibration of the Richards equation with HYDRUS on soil data from the Sunjia catchment. For the six generic soil texture profiles, they are the literature values provided in HYDRUS by the ROSETTA pedotransfer model [29]. For all soil profiles, dtdry, dtsat a, θe and c are the SSMF parameters. The ranges that were used for the calibration of the SSMF are given in the table.
θ r   ( m 3 / m 3 ) θ s   ( m 3 / m 3 ) n (-) l (-) K s   ( m / d ) α (1/m) d t d r y   ( hr ) d t s a t   ( hr ) a (-) c (-) θ e / θ s (-)
CZO clay loam 0–0.3 m0.0500.3901.2560.50.0350.05[0;24][0; d t d r y ][0;100][0;10][0,1]
CZO clay loam 0.3–3 m0.0500.3851.4380.50.0290.086[0;24][0; d t d r y ][0;100][0;10][0,1]
Loam 0–3 m0.0780.4301.560.50.2503.6[0;24][0; d t d r y ][0;100][0;10][0,1]
Loamy sand 0–3 m0.0350.4371.50.51.4664.85[0;24][0; d t d r y ][0;100][0;10][0,1]
Sandy clay loam 0–3 m0.10.391.480.50.3145.9[0;24][0; d t d r y ][0;100][0;10][0,1]
Sandy loam 0–3 m0.0410.4531.3780.50.6213[0;24][0; d t d r y ][0;100][0;10][0,1]
Silt loam 0–3 m0.0670.451.410.50.1082[0;24][0; d t d r y ][0;100][0;10][0,1]
Silty clay loam 0–3 m0.0890.4301.230.51.4661[0;24][0; d t d r y ][0;100][0;10][0,1]

Share and Cite

MDPI and ACS Style

Verrot, L.; Geris, J.; Gao, L.; Peng, X.; Oyesiku-Blakemore, J.; Smith, J.U.; Hodson, M.E.; Zhang, G.; Hallett, P.D. A Simple Modelling Framework for Shallow Subsurface Water Storage and Flow. Water 2019, 11, 1725. https://doi.org/10.3390/w11081725

AMA Style

Verrot L, Geris J, Gao L, Peng X, Oyesiku-Blakemore J, Smith JU, Hodson ME, Zhang G, Hallett PD. A Simple Modelling Framework for Shallow Subsurface Water Storage and Flow. Water. 2019; 11(8):1725. https://doi.org/10.3390/w11081725

Chicago/Turabian Style

Verrot, Lucile, Josie Geris, Lei Gao, Xinhua Peng, Joseph Oyesiku-Blakemore, Jo U. Smith, Mark E. Hodson, Ganlin Zhang, and Paul D. Hallett. 2019. "A Simple Modelling Framework for Shallow Subsurface Water Storage and Flow" Water 11, no. 8: 1725. https://doi.org/10.3390/w11081725

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop