Modelling Soil Water Dynamics from Soil Hydraulic Parameters Estimated by an Alternative Method in a Tropical Experimental Basin

Abstract

Knowledge about soil moisture dynamics and their relation with rainfall, evapotranspiration, and soil physical properties is fundamental for understanding the hydrological processes in a region. Given the difficulties of measurement and the scarcity of surface soil moisture data in some places such as Northeast Brazil, modelling has become a robust tool to overcome such limitations. This study investigated the dynamics of soil water content in two plots in the Gameleira Experimental River Basin, Northeast Brazil. For this, Time Domain Reflectometry (TDR) probes and Hydrus-1D for modelling one-dimensional flow were used in two stages: with hydraulic parameters estimated with the Beerkan Estimation of Soil Transfer Parameters (BEST) method and optimized by inverse modelling. The results showed that the soil water content in the plots is strongly influenced by rainfall, with the greatest variability in the dry–wet–dry transition periods. The modelling results were considered satisfactory with the data estimated by the BEST method (Root Mean Square Errors, RMSE = 0.023 and 0.022 and coefficients of determination, R² = 0.72 and 0.81) and after the optimization (RMSE = 0.012 and 0.020 and R² = 0.83 and 0.72). The performance analysis of the simulations provided strong indications of the efficiency of parameters estimated by BEST to predict the soil moisture variability in the studied river basin without the need for calibration or complex numerical approaches.

1. Introduction

Surface soil moisture plays a key role in the hydrological cycle as it controls the water fluxes between soil, vegetation, and atmosphere [1,2,3]. Moreover, knowledge about soil moisture is widely required in many agricultural studies and applications related to irrigation management [1,4,5,6]. Therefore, monitoring and understanding of soil moisture variability and its exchange relationships with the surface and atmosphere are essential to improve weather forecasting, flooding and drought predictions, and climate projections [7,8]. However, soil moisture is highly variable in space and time due to the combined influence of many factors in a nonlinear fashion such as the hydraulic properties of the soil, topographic characteristics, interaction with surface water systems, precipitation features, and additional meteorological conditions [3,9,10]. Because of its high variability in natural conditions, achieving accuracy in soil moisture estimation to obtain areal information at different spatial and temporal scales is still a challenging task, especially at high spatiotemporal scales [6,11]. In the last 40 years, the scientific community has clearly recognized the importance of soil moisture as input for earth science applications, developing new approaches and techniques for monitoring, modelling, and use of soil moisture data [12].

Currently, three main approaches are used to provide soil moisture estimates: in situ (generally point-scale) measurements, remote sensing observations, and hydrological modelling applications [12,13]. The most accurate and direct methods to determine soil moisture are related to experimental techniques using classical devices such as soil sampling and Time Domain Reflectometry (TDR) probes [14]. However, these measurements are expensive and time-consuming and provide point-scale information but are unable to detect the spatial heterogeneity of the soil moisture at wider scales [6]. The use of remote sensing and modelling techniques has become a viable solution to overcome the limitations of in situ data [15,16]. For instance, remote sensing techniques provide surface soil moisture estimates over large areas by using active sensors operating in microwave bands that are not influenced by solar radiation and cloud cover [6]. However, the spatial and temporal resolutions of the current remote sensing products are still insufficient for small-scale applications (<1 km²) [17], with soil moisture data only available for the first few centimetres of the soil layer [14,18]. Consequently, they do not cover the soil water in the root zone or over the whole profile, which is most interesting for hydrologists [13]. The limitations imposed by the complex interactions between vegetation, soil, and climate make the use of models based on physical processes an effective and viable (low-cost) alternative for evaluating the spatiotemporal soil moisture trends at different depths [2,19,20].

Many models with different levels of complexity and accuracy have been proposed to simulate the dynamics of water in the unsaturated zone of soil; they include, among others, the Soil Water Atmosphere Plant (SWAP) [21], Soil-Vegetation-Atmosphere Transfer (SVAT) [22], Soil Water Infiltration and Movement (SWIM) [23], Groundwater Loading Effects of Agricultural Management Systems (GLEAMS) [24], Water and Agrochemicals in the soil, crop and Vadose Environment (WAVE) [25], and Hydrus models. Most of the abovementioned models use the Richards’ equation to represent the movement of water in unsaturated soil [26], allowing a detailed description of the soil water content distribution and fluxes inside the soil domain [27]. Hydrus-1D model is one of the most common numerical solutions based on the Richards’ equation [28]. Over the years, Hydrus-1D model has been successfully applied in various studies worldwide for predicting soil moisture content and water movement under different conditions (e.g., [13,29,30,31,32,33]). As a physically based model, Hydrus-1D requires some input data for simulations such as meteorological data for surface boundary conditions, soil physical conditions, and physical parameters (saturated hydraulic conductivity and soil water retention curve), which may be obtained through experimentation [27,34,35,36]. Moreover, according to Bordoni et al. [2], for a better implementation of this type of model, a preliminary calibration of some hydrological parameters used in the methodology is sometimes required.

When performed in the laboratory, the methods commonly used to determine soil physical parameters may not be accurate enough to represent the actual field conditions or may be tedious and time-consuming [36,37]. As an alternative, the scientific community has developed simplified methods that are very often able to provide more accurate results, such as Beerkan [38,39]. Beerkan is a simple and easy method based on in situ single-ring water infiltration experiments which can be carried out in the field at low cost [40]. BEST (Beerkan Estimation of Soil Transfer Parameters) is an algorithm presented by Lassabatere et al. [41] which makes it possible to process the infiltration tests carried out using the Beerkan method [42]. The Hydrus-1D model is also able to estimate the soil physical parameters in a simple way through the inverse method, which has been widely used by researchers (e.g., [37,43,44,45]), requiring the use of other types of data observed in the field. However, there is a lack of studies on the performance of Hydrus-1D to simulate the temporal variability of soil moisture from soil physical parameters estimated by the BEST algorithm.

This study focuses on the use of the combination of Hydrus-1D model with BEST method to estimate the temporal variability of soil moisture in an important experimental river basin located near to a highly urbanized area in Northeast Brazil. Specifically, this study aims to: (i) quantify the temporal variability of soil moisture (storage and movement) by using rainfall as the only water input (without irrigation) of the system over two years (2015–2016); (ii) analyze the response of the Hydrus-1D model to simulate the soil moisture dynamics by implementing soil hydrodynamic parameters estimated in the field by the BEST method; and (iii) verify the performance of the Hydrus-1D model after the calibration (inverse modelling) of the soil hydrodynamic parameters. The comparison of modelled data obtained with BEST method and those optimized using numerical inversion allows the characterization of the consistency of BEST estimates for the modelling of water flow. To the authors’ knowledge, such a validation of the use of BEST method has never been done before. Moreover, information about soil moisture dynamics in tropical regions where the monitoring network is scarce are relevant, especially the knowledge of water transfer processes.

2. Materials and Methods

2.1. Study Site Description

This study was carried out in Gameleira Experimental River Basin (GERB), Pernambuco state, Northeast Brazil (Figure 1). The experimental basin covers 17 km² between the coordinates 8°04′ and 8°06′ S and 35°17′ and 35°20′ W. GERB is a sub-basin of the Tapacurá Representative River Basin (TRRB), which plays an important role in water supply and flood control for the Recife Metropolitan Region (RMR) [46], the fourth largest urban area in Brazil, hosting four million inhabitants. The experimental and representative catchments have been monitored by the Network of Hydrology of the Semi-Arid Region (REHISA), which currently consists of eight federal universities in Brazil [47]. In the last years, the hydrology of the GERB has been studied with regard to surface runoff [48], sediment yield [49,50], hydrodynamic soil characterization [51], and evapotranspiration [52,53].

GERB is strategically located near the transition zone between the humid (Atlantic Forest, tropical woodland) and the semi-arid zone (Caatinga, scrub vegetation). The study site is mostly covered by irrigated sugarcane and vegetable crops, evidencing the necessity of practical and precise studies considering information on the water-soil system in the region. Secondary natural vegetation and grassland for extensive cattle-raising are also found in GERB. According to Köppen’s classification proposed by Alvares et al. [54] for the Brazilian territory, the climate in GERB is tropical with dry summer (As), with an average temperature above 23 °C. The average annual rainfall in the experimental basin is approximately 1050 mm, most of which is concentrated from March to July, when the relative air humidity is higher than 70% [46].

The catchment altitude ranges between 140 and 430 m above sea level, with many points presenting high slopes. The dominant soil types in GERB are Acrisols and Gleysols, according to the WRB 2006 [55]. Acrisols are well drained, with structure developed in blocks and a large percentage of silt in their composition. On the other hand, Gleysols occur at the banks of the watercourses and are characerized by a low coefficient of infiltration when saturated and low hydraulic conductivity [56].

2.2. Data Collection

Soil moisture data were automatically measured at two monitoring stations (S1 and S2, Figure 1) by using TDR probes (CS616, Campbell Scientific). Data from monitoring stations S1 and S2 were collected between March 2014 and December 2016. One TDR probe per site was vertically inserted in the soil to indirectly represent the average water content along the first 30 cm downward from the surface, i.e., over the depth of 0–30 cm. The two TDR probes were calibrated in the laboratory, with the reproduction of the field conditions in soil column tests.

Table 1. Physical properties of the soil at the monitoring stations.

Site	Soil Type	Depth (cm)	Soil Particle Composition			Bulk Density	Particle Density
			Sand (%)	Silt (%)	Clay (%)	(g cm−3)
S1	Acrisol	0–10	71.61	20.06	8.33	1.54	2.62
		10–20	71.72	18.68	9.60	1.50	2.60
		20–30	76.14	14.11	9.75	1.51	2.58
S2	Acrisol	0–10	67.56	24.01	8.43	1.58	2.58
		10–20	60.40	21.26	18.34	1.62	2.56
		20–30	65.88	23.10	11.02	1.64	2.58

shows some soil characteristics of monitoring stations S1 and S2, such as the soil type, textural characterization, bulk density, and particle density. Overall, these two sites present similar physical properties of the soil and homogeneity of the profiles along the depths of 0–10, 10–20, and 20–30 cm. During the rainy season, the groundwater levels are approximately 1 and 4 m above the soil surface in stations S1 and S2, respectively. Both monitored sites are covered by sparse grass frequently trimmed, with root depth inferior to 10 cm.

Daily meteorological data were also used in this study to calculate the potential evapotranspiration (ET_p) by Penman’s formulation [57] modified by Shuttleworth [58]. This formulation is a physically-based form of ET_p, which means that all key variables that govern the evaporative process are explicit in the equation [59]. Donohue et al. [60] demonstrated that this form of ET_p is the most appropriated when considering a changing climate, which was also used by Coelho et al. [61] in a study carried out in the Brazilian Northeast region. The data necessary to calculate the ET_p by the Penman’s formulation were acquired from a meteorological station (named S3) operated by Pernambuco Water and Climate Agency (APAC). Meteorological station S3 is approximately 6 km from stations S1 and S2 and automatically monitors the following variables: rainfall, air temperature, relative humidity, atmospheric pressure, wind speed, and solar radiation. In addition, rainfall data were also obtained by using two automatic rain gauges installed in stations S1 and S2.

2.3. Modelling the Unsaturated Water Flow

Unsaturated water flow in GERB was modelled using the Hydrus-1D package, a free modelling environment that allows the simulation of water, heat, and solute transport in variably saturated media under permanent or transient regimes [28]. The governing equation of one-dimensional water flow for a partially saturated porous medium can be described by a modified form of the Richards equation [26], as described in Equation (1). The application of this modified form assumes that the air phase plays an insignificant role in the liquid flow process and the water flow due to thermal gradients, and therefore they can be neglected [62].

$$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left[K\left(h\right)\left(\frac{\partial h}{\partial z}+1\right)\right]-S\left(h\right)$$

(1)

where h is the water pressure head (cm), θ is the volumetric water content (cm³ cm⁻³), t is time (day), x is the vertical coordinate axis (cm) (positive upward), K(h) is the hydraulic conductivity (cm day⁻¹), and S(h) is the sink term in the flow equation (cm³ cm⁻³ day⁻¹), which represents the root water uptake (actual transpiration), obtained using the Feddes equation (more information can be found in the manual [28]). K(h) and θ(h) are obtained from the following functions according to van Genuchten [63]:

$$\theta \left(h\right)={\theta }_r+\frac{{\theta }_s-{\theta }_r}{{\left[1+{\left|\alpha h\right|}^n\right]}^m}, h<0$$

(2)

$$\theta \left(h\right)={\theta }_s, h\ge 0$$

(3)

$$K\left(h\right)=K_s{S}_e^l{\left[1-{\left(1-S_e^{1/m}\right)}^m\right]}^2$$

(4)

where θ_r is the residual water content (cm³ cm⁻³), θ_s is the saturated water content (cm³ cm⁻³), K_s is the saturated hydraulic conductivity (cm day⁻¹), l is the shape factor in the hydraulic conductivity function (assumed to be equal to 0.5), m, α, and n are empirical shape factors in the water retention function (m = 1 − 1/n), and S_e is the relative saturation, calculated as follows:

$$S_e=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}$$

(5)

2.4. Estimation of Soil Hydraulic Parameters Using BEST Methods

The soil hydraulic parameters required for the Hydrus-1D simulation were obtained by the BEST method [64]. The BEST method has wide advantages when compared to other classical methods because it provides a complete characterization of both soil water retention and hydraulic conductivity functions. BEST is based on the van Genuchten [63] relationship for the water retention curve, by using the Burdine [65] condition and the Brooks and Corey [66] relationship for hydraulic conductivity (Equation (6) [65]) and Equation (7). These relationships were used because Fuentes et al. [67] found they gave more accuracy when describing the hydraulic behaviour of most soil types analyzed.

$$\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}={\left[1+{\left(\frac{h}{h_g}\right)}^n\right]}^{-m} \mathrm{with} m=1-\frac{2}{n}$$

(6)

$$\frac{K\left(\theta \right)}{K_s}={\left(\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}\right)}^{\eta }$$

(7)

where h_g is the scale parameter (cm) of θ(h), and η is the shape parameter of the K(θ) relationship. Usually, θ_r is very low and is assumed to be zero in BEST. The saturated water content (θ_s) is derived from the value of bulk density, being assumed to be equal to the porosity. This set of formulae for water retention and hydraulic conductivity functions is referred to as the vGBC functions (van Genucthen–Brooks and Corey functions).

BEST methods require two sets of data to estimate all the hydraulic parameters: (i) the Particle-Size Distribution (PSD) of the soil and bulk density and (ii) cumulative infiltration along with the initial and final soil water contents of the infiltration experiment [40,41]. Three algorithms were developed to analyze the infiltration data: (i) the original version, namely BEST-slope [41], (ii) an improved method for coarse media, namely BEST-intercept [68], and (iii) a third method, BEST-steady [69], based on the analysis of the last part of the cumulative infiltration [40]. The difference between the three methods lies in the way in which the cumulative infiltration data are processed. In this study the BEST-slope and BEST-intercept algorithms were used; however, the results of the first algorithm were discarded due to divergent values of K_s. According to Yilmaz et al. [68], the BEST-slope may lead to erroneous values of K_s, especially when a very high level of precision relative to the steady-state infiltration rate cannot be obtained.

In this study, the infiltration experiments were performed at depths of approximately 10 cm from the soil surface in order to represent the range of soil moisture obtained by the probes installed in the field. For the experiments, fixed water volumes (70 mL) were poured into a 75-mm-diameter ring, leading to a ponded water thickness of a few millimetres. Three infiltration experiments were performed at each of the monitoring stations S1 and S2 to improve the representativeness of the site. The time required for the infiltration of each water volume was measured until steady state was reached and the surface soil was sampled at the end of each experiment. These collected samples were used to determine the final profiles of θ_s and to analyze the particle size of the soils, as previously done by Coutinho et al. [40]. Although BEST assumes that θ_r is equal to zero, the value of this parameter was determined in the laboratory and fixed for the initial model simulations.

The shape parameters n, m, and η are derived from the pedotransfer functions detailed in BEST [41]. Then, the cumulative infiltrations are treated and fitted to the models developed for water infiltration into the disc source by Haverkamp et al. [38] to derive the field-saturated soil hydraulic conductivity (K_s) and soil sorptivity (S). The scale parameter for the water retention curve, hg, is then obtained from the estimates of K_s and S. In this step, all unsaturated hydraulic parameters, namely θ_r, θ_s, n (and thus m), η, h_g, and K_s, are fully estimated and the complete water retention and hydraulic conductivity functions are determined [40,41].

The retention curve parameters obtained by the BEST method were adjusted in this study because the Hydrus-1D model uses the van Genuchten equation with the Mualem [70] pore distribution condition. Therefore, the Retension Curve (RETC) program elaborated by van Genuchten et al. [71] was used. RETC uses the parametric models of vGBC functions to fit the hydraulic curves previously determined with BEST. After adjusting the retention curve, the parameters α = 1/h_g (cm⁻¹) and a new value of n were obtained.

2.5. Boundary Condition

For the atmosphere boundary condition, the variable flow of daily precipitation and potential evapotranspiration were used, considering free drainage for the lower boundary condition. Hydrus-1D requires separated input values of both potential evaporation (the upper boundary condition) and potential transpiration (the sink term in the Richards equation) fluxes at a daily time step to simulate the influence of soil water on transpiration (root water uptake) [62]. Then, potential evaporation (E_p) and potential transpiration (T_p) were separated according to the model based on the Leaf Area Index (LAI) using Beer’s law [72] as follows:

$${\mathrm{T}}_{\mathrm{p}}={\mathrm{ET}}_{\mathrm{p}}\left[1-e^{\left(-k \mathrm{LAI}\right)}\right]$$

(8)

$${\mathrm{E}}_{\mathrm{p}}={\mathrm{ET}}_{\mathrm{p}}{e}^{\left(-0.463 \mathrm{LAI}\right)}$$

(9)

where k is an extension coefficient for global solar radiation and depends on the angle of the sun, the distribution of plants, and the arrangement of leaves [62,73]. In this study, we adopted the same k-value of 0.463 that was considered by Chen et al. [13] and Šimůnek [28]. The value of LAI was assumed to be equal to 1.0 and constant for the whole period of simulation because S1 and S2 are installed in a sparse lawn area. This constant LAI with time is the same used by Chen et al. [13], which selected values in a feasible range for grasses and found LAI = 1 the more consistent after optimization and calibration of the parameters. The partitioned daily values of ET_p were then used as input data for the model (Figure 2).

2.6. Time and Space Discretisation

Regarding the geometry of the numerical domain, we considered a profile with a depth of 60 cm composed of a single control layer (corresponding to the physical and hydrodynamic properties obtained in the infiltration test performed for the 10 cm layer). This assumption of a homogeneous soil along the first 60 cm of the profile was considered after the collection of samples for soil characterization, where the same structural class was found. The numerical domain is made of a compact mesh of 201 nodes. Daily time discretisation was used in this study, with the initial time considered as 0 (the day before the simulation) and the final time referring to the days for calibration and validation. The soil moisture measured on the day before the simulation start was used as the initial condition. For the convergence criteria, the maximum number of 10 iterations for the resolution of the Richards non-linear equation was chosen, with a value of 0.001 for the tolerance of the water content in the unsaturated region of the nodes, in agreement with the default values [28]. As mentioned above, the hydraulic curves were described using the van Genuchten model for the water retention curve [63], with the Mualem condition, and the Mualem model for hydraulic conductivity. This set of formulae for water retention and hydraulic conductivity functions is referred to as the vGM functions. We did not consider the hysteresis phenomenon.

2.7. Simulations of Volumetric Water Content and Model Evaluation

The simulations were started with the hydraulic parameters obtained by the BEST method and adjusted to the vGM formulations for the period of two years. Then, optimization of the hydraulic parameters was performed by inverse modelling with the model parameters adjusted to minimize the difference between observed and simulated values [74]. The parameter algorithm employed by the Hydrus-1D model for optimization purposes is based on the Marquardt-Levenberg method [28]. The optimization was performed by using 365 and 268 days throughout 2015 for S2 and S1, respectively. A reduced number of days was used for this purpose at S1 because the station only started to operate at the end of March 2015.

On the other hand, the year 2016 was used for the validation of the model, performed soon after the calibration of the hydraulic parameters. Three statistical criteria were adopted to analyze the performance of the model: (i) the root mean square error (RMSE); (ii) the Willmott index (d); and (iii) the Nash-Sutcliffe model efficiency coefficient (NSE).

$$\mathrm{RMSE} =\sqrt{\frac{1}{\mathrm{n}}\left[{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\theta }_o-{\theta }_m\right)}^2\right]}$$

(10)

$$\mathrm{d}=1-\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{({\theta }_m-{\theta }_o)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left(\left|{\theta }_m-\overline{\theta }\right|+\left|{\theta }_o-\overline{\theta }\right|\right)}$$

(11)

$$\mathrm{NSE}=1-\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\theta }_o-{\theta }_m\right)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\theta }_o-\overline{\theta }\right)}^2}$$

(12)

where θ_o and θ_m are the observed and modelled soil moisture values, respectively, $\overline{\theta }$ is the mean of the observed moisture values, and n is the number of comparable paired points.

3. Results and Discussion

3.1. Rainfall and Soil Moisture Dynamics

Figure 3 shows the differences between monthly rainfall regimes in the two studied hydrological years and the average monthly rainfall time series from 1970 to 2000 in the Tapacurá representative river basin, TRRB. In 2015, the rainy period occurred between March and July as well as in December. On the other hand, the rainfall was more concentrated from January to May in 2016, which does not correspond to the rainy period registered in the monthly long-term time-series for TRRB [75]. In the two analyzed years, the precipitation in GERB was around 20% less than the annual average for the region, with values of 1821 and 2060 mm registered for S1 and S2, respectively.

The time series of volumetric water content (Figure 4) were obtained based on the calibration curves of CS616 probes performed in the laboratory. Analysing the soil moisture behaviour as a function of precipitation, it is observed that rainfall events above 15 mm generally provide rapid responses in the two stations. The highest increases in soil water content occurred during the rainfall events observed from May to July 2015, with maximum values of 0.32 and 0.33 cm³ cm⁻³ in S1 and S2, respectively. After the abrupt increase of soil moisture at the end of May 2015, which was influenced by a rainfall event of 61 mm, the water content in the soil remained slightly stable in S1 and S2 until the end of the rainy period in July. In 2016, the soil moisture variability was higher than in 2015 during the rainfall events, exhibiting many abrupt increases and decreases in soil water content. These different soil moisture temporal patterns are probably associated with the divergent behaviour of rainfall between the two analyzed years. While in 2015, the rainy season was more concentrated and without the presence of abrupt decreases of soil moisture, throughout the year 2016, the rain was sparser during the rainy season and consequently depletions in the soil moisture were smoother in the two stations. After the rainy period, the water content in the soil presented a mild depletion during the first days of August in both analyzed years, that is, at the beginning of the dry season. These decreases after the rainy season are more pronounced at S2, probably due to the surrounding vegetation slightly lower, and tend to stabilize at values close to 0.20 cm³ cm⁻³.

Figure 5 shows box plots with the variability and distribution of the observed mean daily soil moisture (the dot inside the box) for each month at the two monitoring stations, and the median (central marks), minimum, and maximum are also presented. The box delimits the lower and upper quartiles and outliers. The presence of outliers was only noted in 2015 due to the higher values of soil moisture, probably caused by isolated rainfall events. The isolated rainfall events observed in 2015 generated more accentuated increases in soil water content, differing from the other events. The presence of these atypical values can also be associated with the beginning and ending of the rainy season, as observed in May and August.

The highest dispersion of data is observed in the transition months between dry–wet and wet–dry periods, such as August and December 2015 and the first semester of 2016, and during the wetter periods. This variability can be noted from the difference between maximum and minimum values as well as the distance between the quartiles and the median. On the other hand, the driest months presented the lowest variability of data, as expected due to the absence of high rainfall volumes. These lowest variabilities were observed from September to October 2015 and after the end of the rainy season in 2016 (starting in August).

3.2. Hydraulic Characterization of the Sites

Table 2. Mean values of the hydraulic parameters derived from Beerkan Estimation of Soil Transfer (BEST) and Retension Curve (RETC) of the three replicates performed at each station.

Site	BEST Output								RETC Output
	m	n	η	S (mm s−1/2)	Ks (cm day−1)	θr1 (cm3 cm−3)	θs (cm3 cm−3)	hg (cm)	A (cm−1)	n
S1	0.102	2.223	11.482	1.31	159.97	0.063	0.377	−18.5	0.041	1.245
S2	0.094	2.204	12.682	1.20	155.52	0.026	0.379	−16.7	0.020	1.273

¹ The value of θ_r was measured in the laboratory.

presents the results of the mean values of the shape (m, n, and η) and normalisation (S, K_s, θ_s, and h_g) parameters of the soils generated by the BEST method for stations S1 and S2. Adjusted parameters (α and n) for the vGM functions using the RETC program are also exhibited in Table 2. The parameters m and n, which depend on the soil texture, presented the highest values in the soil with the sandiest fraction, that is, at station S1. The shape parameters strongly depend on the type of soil, with higher values of parameters n and m for the coarser materials (sand and silty sand) and smaller values of the parameter η for coarse soils [40,76]. The soil at station S1 also exhibits slightly higher values of S and K_s. The estimated values are similar to the results found by Furtunato et al. [51], who also used the BEST method at 102 points distributed throughout the GERB.

Figure 6 shows the hydraulic conductivity and water retention curves of the soils at stations S1 and S2 obtained from the mean values of the hydraulic parameters. The water retention curves were fitted to the vGM functions, generating the new parameters α and n (Table 2).

Accordingly, a similar hydraulic behaviour was verified between the soils present at the studied stations. Both hydraulic curves present an accentuated format typical of coarser materials, with a very marked inflexion point before saturation. This similarity of shapes is in agreement with the proximity of the adjusted values of the parameter n in the two soils (Table 2), as this parameter mainly determines the slope of the retention curve at larger potential values [77]. In this case, water retention is mainly affected by the soil texture [78], which is similar in the two stations (see Table 1). It is possible to observe that the Ks values at station S1 are higher than those noted at S2 for almost all ranges of soil water content. The values of K at the two stations only reach similar amounts when the soils are near to saturation (~ 160 cm day⁻¹).

3.3. Simulation

The first simulations were performed with the hydraulic parameters estimated by the BEST method and adjusted to vGM functions (Table 2). The results of the simulation at both station S1 and station S2 present similar tendencies of variation, with estimated soil moisture values close to the observed data in almost the entire series throughout the studied period (Figure 7). Exceptions are noted during the dry periods, especially in S2, when the simulated data cannot effectively capture the observed recession curves. This may have been influenced by the lower boundary condition adopted in the simulation, which in this case was free drainage disregarding the depth of the groundwater. It is possible that in the field there is a small contribution of capillarity to soil moisture. Another possible explanation for the underestimated results of the model at the two sites may be the constant value considered for LAI over the whole simulation, promoting periods with higher evaporation rates and consequently reducing soil moisture rates.

The scatter plots also show the good agreement between observed and estimated data, with values of R² ranging from 0.72 (S1) to 0.81 (S2) and RMSE below 0.024 cm³ cm⁻³ (Figure 7b,d). Honari et al. [36] used the soil parameters estimated by the BEST method to evaluate the performance of Hydrus-3D model in simulating 3D water flow and water content in a short time period and at different depths of a soil cultivated with corn and durum in Montpellier, France. They also obtained good results without parameter calibration, with RMSE values between 0.0174 and 0.019 cm³ cm⁻³. The authors emphasized that the quality of hydrodynamic property estimates is fundamental for accuracy in simulations.

Some studies used hydraulic parameters estimated by the Rosetta Lite pedotransference function [79] to start the simulations of the behaviour of water in the vadose zone using Hydrus-1D [30,32,33,73]. By analysing the efficiency of this function, Okamoto et al. [80] observed that the hydraulic parameters measured in the field presented better results than those estimated by Rosetta Lite. They highlighted the importance of the use of measured parameters when simulating the water flux in the unsaturated zone. The results found at this step demonstrate the good efficiency of hydrodynamic parameters determined in situ with the BEST model, adjusted by the RETC, to characterize the temporal variability of soil moisture in the GERB.

Table 3. Calibrated soil hydraulic parameters by using the inverse solution for the monitoring stations.

Site	θr (cm3 cm−3)	θs (cm3 cm−3)	α (cm−1)	n	Ks (cm day−1)	l
S1	0.107	0.342	0.037	1.23	234.2	0.5
S2	0.056	0.369	0.018	1.17	243.7	0.5

presents the hydraulic parameters of the soil calibrated by inverse modelling. After calibration, residual moisture values increased and saturated moisture decreased in the two stations. The saturated hydraulic conductivity in the two stations increased approximately 50% of predicted values using the BEST methods (Table 2 versus Table 3). Chen et al. [13] used Hydrus-1D to simulate soil moisture in two sub-basins of the Goulburn River, Australia, for a period of three years. They observed that inverse modelling modified the soil properties, reducing the values of θ_r while θ_s and K_s increased. According to the aforementioned authors, the calibrated hydraulic parameters represent, in addition to soil texture, factors such as the number of macropores caused by roots or cracks, vegetation type, organic matter, and anthropic activities.

Figure 8 shows the results of the modelling validation. It is noticeable that the estimated soil water content is better adjusted to the observed values obtained in situ. However, a discrepancy between the observed and estimated data remains in both stations after the end of the rainy season. At S1, the simulated data follow the same pattern of the measured soil moisture variation; however, the model slightly overestimated the measured values for some rainfall events (Figure 8a). For instance, the response of the model to the rainfall of 133.8 mm that occurred on 9 May 2016 was 0.33 cm³ cm⁻³, while the observed value registered was 0.27 cm³ cm⁻³. For the same station, the model was unable to detect the rainfall events that occurred from 13 June 2016 to 10 August 2016, underestimating the measured soil moisture data.

The simulated data for S2 also presented overestimated values for some rainfall events (Figure 8). Similar underestimation was also observed after the rainy season from 10 July 2016 to 9 August 2016, with a considerable difference between the observed and estimated data. This behaviour of the model in underestimating the observed data is probably linked to the high infiltration rate considered by Hydrus-1D, which is larger than that which actually occurs with the antecedent soil moisture.

The correlation between estimated and simulated values for the two monitoring stations is shown in Figure 8b,d. The dispersion between the data indicates that the observed soil water content is represented well by the modelling results, with R² values equal to 0.83 (S1) and 0.72 (S2). Silva et al. [29] obtained similar values of R² (0.79 for natural cover and 0.78 for bare soil) by using the Hydrus-1D model to simulate soil moisture in an Acrisol in Brazil using 151 days of the year 2013 for validation. The authors attributed the results obtained to the different hydrological behaviours between the years used for calibration and validation, particularly the changes in interception rate and evapotranspiration and the scarcity of rainfall events in the calibration period. Any change in vegetation or soil structure at the surface may impact the water budget, mainly the runoff and infiltration processes.

Figure 8b,d also shows that some points are still far from the straight 1:1 line for the two stations, even though R² presents good values. This occurred because of the response of the model to some isolated rainfall events, which raised the soil moisture to values higher than those measured in the field. Another divergent behaviour is observed below the 1:1 line, where the model failed to respond to the small rainfall events that occurred in June and July after the rainy season, underestimating the soil moisture values measured in the field.

Another possible cause of the observed divergence between simulated and observed values is associated with some possible small errors that accumulated throughout the data acquisition process, such as: (i) characterization of soil properties, (ii) calibration of the sensors in the laboratory, and (iii) some criteria adopted for simulation of the model (e.g., LAI and lower boundary condition), as previously described.

Table 4. Hydrus-1D model goodness of fit and accuracy. RMSE: Residual Mean Square Errors; d: Willmott coefficient of agreement; NSE: Nash-Sutcliffe model efficiency coefficient.

Site	RMSE (cm3 cm−3)	d	NSE
S1	0.012	0.95	0.783
S2	0.020	0.93	0.740

presents the results of the performance evaluation of the model. The RMSE for the simulated and observed soil moisture data exhibited good absolute values equal to 0.012 and 0.020 cm³ cm⁻³ for S1 and S2, respectively. The Willmott coefficient of agreement (d) close to 1 also evidences the good accuracy of the soil moisture simulations under the investigated conditions. The NSE values for the simulations at both station S1 and station S2 were classified as being “very good” (0.75 < NSE < 1) and “good” (0.65 < NSE ≤ 0.75), respectively. These results are similar to those found by Silva et al. [29] using Hydrus-1D for simulations, with d equal to 0.93 and 0.94 and NSE equal to 0.72 and 0.74 for natural and bare soils, respectively.

The comparison between simulated and validated results shows that the RMSE values were very close, with a small reduction of 0.009 cm³ cm⁻³ at S1 and a of 0.002 cm³ cm⁻³ at S2 (see Figure 7 and Figure 8). Considering the aforementioned results and the entire time series used in both the initial simulation and the validation (2015 and 2016), it can be noticed that the parameters estimated by the BEST method are efficient for use in the Hydrus-1D model and describe the soil conditions in the GERB well, without the need for calibration. Therefore, simplified estimates of soil properties by the BEST method, associated with meteorological data and using Hydrus-1D model, can provide acceptable predictions of soil water content in GERB.

4. Conclusions

This study was conducted in an experimental basin in Brazil to investigate the behaviour and prediction of soil water content dynamics in the region using field measurements and numerical methods. Soil moisture presented the greatest and the lowest variability during the rainy transitions (dry–rainy–dry) and dry periods, respectively. The soil water content did not approximate the residual and saturation contents throughout the period. The Hydrus-1D model adequately simulated the soil moisture variability both with the parameters estimated by the BEST model and after their calibration (numerical optimization of parameters by fitting experimental measures). Indeed, these optimized parameters reduced the differences between the measured and simulated values, as expected. However, the inverse modelling requires a considerable number of observations, which means it is limited to well monitored basins. Meanwhile, the gain in accuracy was not drastic, suggesting that BEST estimates were sufficient. Thus, the performance analysis of the simulations provided strong indications of the efficiency of parameters estimated by BEST to predict the soil moisture variability in the studied region, without the need for calibration or complex numerical approaches. This result indicates that the BEST method can be adequate for depicting the water retention and hydraulic conductivity functions, even with the final aim of modelling the water flow processes and water budget for a soil profile. It must be mentioned that the approach was based on the modelling of a soil profile, considering a uniform soil. The impact of soil heterogeneity and spatial variability of soil hydraulic characteristics in different layers of the soil profile should be investigated and the resulting flow heterogeneity should be accounted for. This issue will be the object of further investigations.