A Simple Modelling Framework for Shallow Subsurface Water Storage and Flow

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.


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 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.

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].
he Sunjia CZO site and 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 USD" 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.

Model Requirements
The SSMF was developed to meet the following requirements:

1.
Its results should be comparable to the Richards equation solution; 2.
It should be computationally simple and parsimonious, with a minimum number of calibrating parameters; 3.
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 (m 3 /m 3 ) as well as soil water flux L d (e.g., m/d) across a shallow depth z r (m). 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, K s 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 z ls (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.

Model Conceptualization and Development
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 (m 3 /m 3 ), above which a negative surface water balance (P − E < 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 dt sat (hr), and longer when the soil is dry, reaching a value of dt dry (hr). Between those two values the time step follows a linear relationship. For a daily time-scale, the minimum and maximum values for dt sat and dt dry vary between 0 and 24 h, respectively, with dt dry > dt sat . The values of dt sat and dt dry 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 ls d,pot (m/d), computed from the daily total evapotranspiration ET d (m/d) and the total daily precipitation P d (m/d) as shown in Equation (1): where ET 0,d is daily total potential evapotranspiration (m/d), θ d−1 (m 3 /m 3 ) is the soil water content of the previous day, as calculated in the following routine (Equation (8)), θ r is the residual water content (m 3 /m 3 ), θ s is the saturated water content (m 3 /m 3 ) 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, E w , (evapotranspiration at the wilting point) from the original publication is omitted here for simplicity. Rather than using this method to obtain ET d , the SSMF used Penman-Monteith to estimate ET 0,d . The daily subroutine was then divided into consecutive time steps t by the following equations. We first determined how much of q ls t,pot infiltrated into the soil at the time step t (q ls t (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)).
This input of water led to a temporary soil water content θ t,1 (m 3 /m 3 ), as defined in Equation (4): θ 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: for which dt sat (hr), and dt dry (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, dt dry must be relatively high, while dt sat must be relatively low, between 0 h and dt dry . 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 t tot (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 z r t,down (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): where θ s (m 3 /m 3 ) is the saturated water content, θ r (m 3 /m 3 ) 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 z r t,down was then defined as 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=t−1 k=0 q z r k,down dk mitigates the vertical flux through z r . 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=t−1 k=0 q z r k,down dk 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 z r t,down . 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 : Finally, the potential soil water recharge from the land surface for the next time step q ls t+1,pot was updated The sub-daily routine stopped when the condition t 0 dt i = t tot was met. The daily depth-averaged soil-water content value θ d and vertical downwards flux across z r , q z r d,down were, therefore, the values at the last time step of the sub-daily routine (t = t end ) and the sum of the downwards flux across z r over the sub-daily routine, respectively, calculated as: The vertical upwards flux across z r q z r d,up (Equations (12a) and (12b)) was triggered when the soil water content was large, above θ e , and when 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: Eventually when the total time t tot was reached, the excess water constituted the overland flow where ∆S d (m/d) is the change in the soil water storage over z r ; For the first day of the modelling time period, d = 1, the soil water content was equal to the initial value 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 ls t=0,pot = q ls d,pot 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 ls t=0,pot = 0. The initial value of the soil water content was the soil water content of the previous day (Equation (16)):

Model Evaluation
For all seven soils, HYDRUS was used to derive daily soil water content θ d , HYDRUS and soil water flux dynamics L d,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,HYDRUS . Next, the SSMF was calibrated against the soil water content θ d , HYDRUS through a Monte-Carlo simulation (10,000 runs) and its results L d,mod were compared with L d,HYDRUS. Figure 2 provides a summary of this procedure and Table 1 provides the input soil parameters θ r , θ s , n, l, K s , α and SSMF parameter calibration ranges for dt dry , dt sat and θ e (as a value relative to the saturated water content θ s , a and c. 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, K s 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, dt dry , dt sat a, θ e and c are the SSMF parameters. The ranges that were used for the calibration of the SSMF are given in the table. 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)).
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 dt dry and dt sat , the SSMF was computed with a constant hourly timestep (dt t = dt = 1 h, Equation (5)).
The potential evapotranspiration ET 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 ET d calculated by the SSMF cannot be directly compared to the value obtained from HYDRUS. This is because ET 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 ET d could be implemented, if those fluxes were the focus of interest.

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.  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 z r = 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 dt dry ; dt sat 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 z r (from NRMSD = 0.09 at z r = 0.10 m to NRMSD = 0.06 at z r = 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 z r regardless of soil, particularly as z r increased.

SSMF Parameters Influence
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 z r . 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).

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 (Figures 4, 5 and 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 z r . 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 z r and that may not represent accurately the properties lower than z r .
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.

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 (ET d ) exceeds the precipitation (P d ), as would occur during dry periods, setting θ e to 0 results in an upwards flux, q z r d,up , that equals the total potential flux (q ls d,pot ). The upwards flux recharging z r 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 dt dry , dt sat , 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 z r (Equation (7)), so the vertical flux equals the hydraulic conductivity corresponding to the conditions in the upper zone z r (unit-hydraulic gradient condition) . This causes the differential in soil water content between the upper zone and the soil below z r to be large and positive, producing wetter conditions in the upper soil than in the soil below z r 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 z r (Equation (7)) for wet water contents.

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 (dt dry and dt sat , a, and θ e ) together lead in almost all the soil profiles and all values of z r 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.

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 http://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.