The Distributed Xin’anjiang Model Incorporating the Analytic Solution of the Storage Capacity Under Unsteady-State Conditions

Abstract

Developing a functional linkage between hydrological variables and easily accessible terrain and soil information is a novel concept for distributed hydrological models. This approach aims to address limitations imposed by data scarcity and high computational demands. The model hypothesizes that the relationship between the evaporation flux and the absolute value of the matric potential follows a power exponential pattern. Analytic solutions for the groundwater depth, the evaporation capacity, and the storage capacity are derived with respect to the topographic index, considering the relationship between the groundwater depth and the topographic index and the influence of setting off. Subsequently, a distributed Xin’anjiang Model using the analytic solution of the storage capacity under unsteady-state conditions is constructed. This new model is employed to simulate soil moisture and discharge in the Tarrawarra Watershed. The simulation results for soil moisture and discharge are compared with those from the Storage Capacity Model and the DHSVM. Additionally, the computational speeds of all three models are compared. The findings indicate that the simulation accuracy of the new model for soil moisture and discharge surpasses that of the Storage Capacity Model and the DHSVM. Meanwhile, the computational speed of the new model is significantly faster than the DHSVM and slightly slower than the Storage Capacity Model. It offers a balance between computational efficiency, predictive accuracy, and physical mechanism representation. The data requirements of the new model are minimal and easy to procure, and it requires less computational effort. Moreover, it accurately captures the spatial and temporal dynamics of soil moisture and the discharge process of the watershed.

1. Introduction

Distributed hydrological models excel in capturing the spatial and temporal dynamics of hydrological variables [1,2,3,4], playing a pivotal role in hydrology and water resource simulation, ecological impact assessment, and hydrological forecasting in regions with limited data [5,6,7,8,9,10]. These models typically fall into two categories: one developed through the numerical solution of microcell hydrophysical equations, exemplified by models such as SHE [11] and DHSVM [12]. However, these models face challenges due to their complex calculations and the need for extensive parameter calibration. The uncertainty in parameter determination, based on limited observational data, poses a significant challenge. Moreover, the computational intensity required for calibrating parameters across numerous grid points can be prohibitive, constraining the application of such models due to the substantial computational resources needed [13].

During the past two decades, the development of distributed hydrological models has garnered significant interest due to their ability to establish functional relationships among hydrological variables and terrain, soil, and other surface information that is readily accessible. Key models in this field include the distributed Xin’anjiang Models, which were proposed by researchers such as Peng Shi [14], Xi Chen [15], Cheng Yao [16], and Xiaohua Xiang [17]. These models have made significant strides in understanding the hydrological processes by incorporating topographic and soil characteristics. The distributed Xin’anjiang Models have been enhanced with analytic solutions for storage capacity related to topographic indices, allowing for the determination of spatial distribution of storage capacity and thus endowing the model with distributed functionality [18]. However, these models have made certain simplifications regarding the hydrodynamic characteristics of the unsaturated zone, which may have impacted the physical mechanisms and introduced some level of distortion in the calculations. Peng Shi [14] proposed an empirical relationship between the storage capacity and the topographic index using a log Weibull function. Xi Chen [15] and Cheng Yao [16] assumed that the soil moisture in the unsaturated zone equates to the wilting moisture content. Xiaohua Xiang [17], on the other hand, disregarded the impact of evaporation on soil water potential in the unsaturated zone, implying a vertical equilibrium state. These simplifications, while making the models more tractable, may have reduced their physical fidelity. Sadeghi, in an attempt to address these limitations, considered the effect of evaporation and derived an analytic solution for water deficit related to matric potential [19]. This approach provided a more nuanced understanding of water deficit calculations. However, Sadeghi’s model still assumed a steady state in the vertical direction and simplified the evaporation flux, which could lead to discrepancies between the calculated and actual soil vertical profiles. In summary, while the distributed Xin’anjiang Models and the work of these researchers have contributed to advancements in hydrological modeling, there is a need for further refinement to better capture the complex dynamics of the unsaturated zone and to reduce the errors in model predictions. This could involve more detailed representations of soil moisture dynamics, improved understanding of the effects of evaporation, and the incorporation of more sophisticated methods for modeling the vertical soil profile.

In this scholarly work, a novel Xin’anjiang Model is introduced, which employs an analytic solution for assessing the storage capacity under unsteady-state conditions. This advancement is a product of extensive research by hydrologists worldwide. The model postulates a power exponential relationship between the evaporation flux and the absolute value of matric potential. The model integrates the impact of setting off by assuming a base storage capacity (W0) at the wettest grid point (referred to as the fundamental unit derived from a grid-type Digital Elevation Model) and designates the maximum groundwater depth at this point as z. By correlating the groundwater depth with a topographic index, analytic solutions for the groundwater depth, the evaporation capacity, and the storage capacity are deduced. The parameters of the power exponential function, denoted as M and x, are derived using the aforementioned analytic solutions, along with the phreatic evaporation. By utilizing the distribution of the topographic index, one can derive the distribution of the storage capacity under unsteady-state conditions. This enables the distributed computation of the new Xin’anjiang Model, which, in theory, offers high accuracy in simulating hydrological variables. Nevertheless, it is essential to validate the computational accuracy and processing speed of the new model introduced in this paper to confirm its practical applicability and reliability in the field of hydrological research.

This study selects discharge data, meteorological data from one weather station, and soil moisture data from 17 observation points within the Tarrawarra Watershed in Australia to calibrate and validate the parameters of a newly developed hydrological model. The objective is to ascertain the computational precision of the new model by evaluating its simulation of spatiotemporal soil moisture variations and the discharge dynamics of the watershed against those of two established models: the Storage Capacity Model (proposed by Xiaohua Xiang as the distributed Xin’anjiang Model) and the DHSVM (Distributed Hydrological Soil-Vegetation Model). Furthermore, the computational efficiency of the new model is benchmarked against the other two models by comparing the time required for 100 stochastic parametric simulations.

2. Materials and Methods

2.1. Study Area

The Tarrawarra Watershed, situated on the periphery of Melbourne, Australia, serves as the focal point of this research. Spanning an area of 10.5 hectares, the watershed is characterized by a mild climate and receives an average annual precipitation of 820 mm [20]. The soil exhibits a distinct seasonal moisture pattern, with wet conditions prevalent in winter and dry conditions in summer. There is the situation of water flow interruption in the watershed during some parts of the year because the area of the watershed is small. The predominant vegetation in the Tarrawarra Watershed consists of grasses. The soil in this region is predominantly composed of silt and clay, with a composition that is approximately 65% silt and 35% clay. This particular blend is commonly designated as silty clay soil. Soil moisture variation occurs within a depth of 0 to 60 cm below the ground surface. Notably, the soil moisture fluctuations in the upper 30 cm layer contribute significantly, accounting for 40% to 60% of the total variation observed in the watershed [21].

The Tarrawarra Watershed data set encompasses a comprehensive range of hydrological and meteorological data, including 5 m contour terrain data (as depicted in Figure 1), discharge data, and meteorological data from a single weather station. Additionally, the data set comprises soil moisture data from 17 observation points and groundwater level data from 74 observation points. Soil moisture is measured at six distinct depth intervals: 15 cm, 30 cm, 45 cm, 60 cm, 90 cm, and 120 cm below the ground surface (using a soil moisture content tester to measure the soil water content, Shanghai Yuanxin Biotechnology Company, Shanghai, China). Observations are recorded at an interval of 8 days. Among the soil moisture data, the highest recorded value is 48.2%, observed at site S9, while the lowest is 21.2%, recorded at site S10. These data points are instrumental in validating the simulation accuracy of the distributed hydrological model, particularly in terms of the temporal and spatial dynamics of hydrological variables. Furthermore, the continuous discharge data serve to verify the model’s simulation effectiveness concerning the watershed’s discharge.

Terrain data were acquired from remote sensing systems, while discharge, soil moisture, and groundwater depth data were collected through surface observations. Meteorological data were obtained from a combination of surface observations, meteorological satellites, radar, and automatic stations. The accuracy and representativeness of these data sets have been previously discussed in References [20,21]. Based on these discussions, it can be concluded that the data utilized in this study are reliable.

The distribution of topographic indices within the Tarrawarra Watershed has been determined utilizing the 5 m contour terrain data. The analysis reveals a minimum topographic index of 4.24, a maximum of 11.91, and an average value of 6.45, as illustrated in Figure 2.

2.2. Data Collection

All data pertinent to the Tarrawarra Watershed, including the 5 m contour terrain data, discharge data, meteorological data from a single weather station, soil moisture data from 17 observation points, and groundwater level data from 74 observation points, have been procured from the designated website: http://people.eng.unimelb.edu.au/aww/tarrawarra/datapage.html (accessed on 1 July 2024)

Furthermore, the topographic indices for the watershed have been calculated based on the 5 m contour terrain data available from the same source.

The terrain data employed in this study originates from 5 m contour maps. Conventional topographic data typically have a resolution of 30 m. However, given the relatively small size of the watershed under investigation, it is crucial to use higher resolution terrain data to ensure the accuracy of the calculations when utilizing this data to determine various variables. The higher resolution allows for a more precise representation of the watershed’s topography, which is essential for accurate modeling and analysis. Therefore, we have chosen to use the 5 m contour terrain data in this paper to meet these requirements.

2.3. Methods

2.3.1. The Analytic Solutions of Groundwater Depth and Evaporation Capacity About Matric Potential Under Unsteady-State Conditions

Numerous models have been developed to describe soil moisture dynamics. Prominent among these are the Brooks–Corey model [22], the Van Genuchten model, and its modified version [23]. Each of these models demonstrates a high degree of accuracy in simulating soil moisture conditions. In the present study, we have chosen to utilize the Brooks–Corey model for the derivation of our mathematical formulas.

Based on the Brooks–Corey model, the variables can be expressed as

$${\displaystyle \frac{K(\psi )}{Ks}}={({\displaystyle \frac{\psi }{{\psi }_b}})}^{-p}, {\displaystyle \frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}}={({\displaystyle \frac{\psi }{{\psi }_b}})}^{-\lambda }, {\displaystyle \frac{{\theta }_f-{\theta }_r}{{\theta }_s-{\theta }_r}}={\left({\displaystyle \frac{{\psi }_c}{{\psi }_b}}\right)}^{-\lambda }$$

(1)

where $K(\psi )$ signifies the hydraulic conductivity coefficient (mm/day); $Ks$ signifies the saturated hydraulic conductivity coefficient (mm/day); $\psi$ denotes the absolute value of the matric potential (mm); ${\psi }_b$ denotes the intake potential (mm); $\theta$ denotes the soil moisture content; ${\theta }_r$ denotes the wilting moisture content; ${\theta }_f$ denotes the field moisture capacity; ${\theta }_s$ denotes the saturated moisture content; ${\psi }_c$ denotes the absolute value of the matric potential, which corresponds to the field moisture capacity; and $p,\lambda$ represent the parameters of Brooks–Corey model, with the relationship $p=2+3\lambda$. Additionally, it is defined such that ${\theta }_f-{\theta }_r=0.85({\theta }_s-{\theta }_r)$.

Reference [24] demonstrates the use of an empirical function to simulate the concentration of Cs in the river, yielding promising results. Drawing from this successful application, we propose to employ a similar empirical function to model the relationship between the evaporation flux and the absolute value of the matric potential.

For a particular grid point characterized by a specific groundwater depth, it is assumed that the relationship between the evaporation flux and the absolute value of the matric potential can be mathematically articulated as

$$E=M{\psi }^x$$

(2)

where $E$ signifies the evaporation flux (mm/day); M, x represent the parameters of the power exponential function that models the evaporation flux as a function of the absolute value of the matric potential; and $\psi$ denotes the absolute value of the matric potential (m).

The diving evaporation can be expressed as

$$E_g=E_p\cdot {(1-{\displaystyle \frac{d}{d_{\max }}})}^{nn}$$

(3)

where $E_g$ signifies the diving evaporation (mm/day); $E_p$ denotes the evaporation capacity (mm/day); $d$ corresponds to the groundwater depth (mm); $d_{\max }$ represents the extreme depth of the diving evaporation (mm); and $nn$ signifies the parameter of the diving evaporation.

According to the Richards equation [25], when $\psi >{\psi }_b$, it can be expressed as

$$Ks{({\displaystyle \frac{\psi }{{\psi }_b}})}^{-p}({\displaystyle \frac{\partial \psi }{\partial z}}-1)=M{\psi }^x$$

(4)

$$\mathrm{Define} {\psi }_m={({\displaystyle \frac{K_s{{\psi }_b}^p}{M}})}^{1/(x+p)}, {\psi }_m{\displaystyle \frac{1}{1+{(\psi /{\psi }_m)}^{p+x}}}d(\psi /{\psi }_m)=dz$$

(5)

$$\begin{gathered} \mathrm{Define} {\displaystyle \frac{\psi }{{\psi }_m}}=y, \mathrm{when} y<1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ {\displaystyle \frac{1}{1+y^{p+x}}}={\displaystyle \sum _{i=0}^N{(-1)}^i}{y}^{i(p+x)}, {\displaystyle \int \frac{1}{1+y^{p+x}}}dy={\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{y}^{1+i(p+x)}}{1+i(p+x)}} \end{gathered}$$

(6)

$$\begin{gathered} \mathrm{When} y>1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ {\displaystyle \frac{y^{-(p+x)}}{y^{-(p+x)}+1}}={\displaystyle \sum _{i=1}^N{(-1)}^{i+1}}{y}^{-i(p+x)}, {\displaystyle \int \frac{y^{-(p+x)}}{y^{-(p+x)}+1}}dy={\displaystyle \sum _{i=1}^N\frac{{(-1)}^{i+1}{y}^{1-i(p+x)}}{1-i(p+x)}} \end{gathered}$$

(7)

$$\begin{gathered} \mathrm{If} {\psi }_m>{\psi }_b\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ z(\psi )=z({\psi }_b)+{\psi }_m{\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{(\psi /{\psi }_m)}^{1+i(p+x)}}{1+i(p+x)}}-{\psi }_m{\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{({\psi }_b/{\psi }_m)}^{1+i(p+x)}}{1+i(p+x)}} \\ {\psi }_b<\psi <={\psi }_m \end{gathered}$$

(8)

$$\begin{gathered} z(\psi )=z({\psi }_m)+{\psi }_m{\displaystyle \sum _{i=1}^N\frac{{(-1)}^{i+1}{(\psi /{\psi }_m)}^{1-i(p+x)}}{1-i(p+x)}}-{\psi }_m{\displaystyle \sum _{i=1}^N\frac{{(-1)}^{i+1}}{1-i(p+x)}} \\ \psi >{\psi }_m \end{gathered}$$

(9)

$$\begin{gathered} \mathrm{If} {\psi }_m<={\psi }_b\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\psi >{\psi }_b\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ z(\psi )=z({\psi }_b)+{\psi }_m{\displaystyle \sum _{i=1}^N\frac{{(-1)}^{i+1}{(\psi /{\psi }_m)}^{1-i(p+x)}}{1-i(p+x)}}-{\psi }_m{\displaystyle \sum _{i=1}^N\frac{{(-1)}^{i+1}{({\psi }_b/{\psi }_m)}^{1-i(p+x)}}{1-i(p+x)}} \end{gathered}$$

(10)

Let $z({\psi }_s)=D$ denote the depth of groundwater. The absolute value of the matric potential, represented by ${\psi }_s$, corresponding to $z({\psi }_s)$, can be determined using a trial algorithm (${\psi }_s$ is also interpreted as the absolute value of the surface matric potential).

Where ${\psi }_s$ signifies the absolute value of the surface matric potential (mm); $z({\psi }_s)$ signifies the vertical distance from the groundwater table to the surface (mm); and $z({\psi }_b)$ denotes the vertical distance from the groundwater table to the point where the absolute value of the matric potential is ${\psi }_b$ (mm).

The relationship between the actual evaporation and the soil moisture content can be described as follows (according to the Xin’anjiang Model [18]):

$$E=E_p \theta >{\theta }_f, E=E_p\cdot (a1+(1-a1){\displaystyle \frac{\theta -{\theta }_r}{{\theta }_f-{\theta }_r}}) {\theta }_r\le \theta \le {\theta }_f$$

(11)

where $E$ signifies the actual evaporation (mm/day) and $a1$ denotes the parameter concerning soil evaporation.

The evaporation capacity can be expressed as

$$\begin{gathered} \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{E}_p={\displaystyle {\int }_{E_c}^{E_s}d{E}_p}+M{\psi }_c^x={\displaystyle {\int }_{{\psi }_c}^{{\psi }_s}\frac{Mx{\psi }^{x-1}}{a1+(1-a1)\frac{\theta -{\theta }_r}{{\theta }_f-{\theta }_r}}d\psi }+M{\psi }_c^x \\ ={\displaystyle \frac{Mx}{a1}}{\displaystyle {\int }_{{\psi }_c}^{{\psi }_s}\frac{{\psi }^{x-1}d\psi }{1+\frac{(1-a1){\theta }_s{\psi }^{-\lambda }}{a1{\theta }_f{\psi }_b^{-\lambda }}}}+M{\psi }_c^x \end{gathered}$$

(12)

where ${\psi }_c$ signifies the absolute value of the matric potential corresponding to the field moisture capacity (mm); $E_c$ denotes the evaporation capacity under the point where the matric potential is equivalent to the field moisture capacity (mm/day); and $E_s$ denotes the evaporation capacity under the surface (mm/day).

Define ${\psi }_g={\left({\displaystyle \frac{a1{\theta }_f}{(1-a1){\theta }_s{\psi }_b^{\lambda }}}\right)}^{-1/\lambda }$, Formula (12) can be expressed as

$$E_p={\displaystyle \frac{Mx{\psi }_g^x}{a1}}\cdot {\displaystyle {\int }_{{\psi }_c}^{{\psi }_s}\frac{{(\psi /{\psi }_g)}^{x-1}d(\psi /{\psi }_g)}{1+{(\psi /{\psi }_g)}^{-\lambda }}}+M{\psi }_c^x$$

(13)

$$\begin{gathered} \mathrm{Define} y={\displaystyle \frac{\psi }{{\psi }_g}}, \mathrm{when} y<1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ {\displaystyle \frac{y^{\lambda +x-1}}{1+y^{\lambda }}}={\displaystyle \sum _{i=0}^N{(-1)}^i}{y}^{i\lambda +\lambda +x-1}, {\displaystyle \int \frac{y^{\lambda +x-1}}{1+y^{\lambda }}}dy={\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{y}^{(i+1)\lambda +x}}{(i+1)\lambda +x}} \end{gathered}$$

(14)

$$\begin{gathered} \mathrm{when} y>1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ {\displaystyle \frac{y^{x-1}}{y^{-\lambda }+1}}={\displaystyle \sum _{i=0}^N{(-1)}^i}{y}^{-i\lambda +x-1}, {\displaystyle \int \frac{y^{x-1}}{y^{-\lambda }+1}}dy={\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{y}^{-i\lambda +x}}{-i\lambda +x}} \end{gathered}$$

(15)

$$\begin{gathered} \mathrm{If} {\psi }_g<={\psi }_c\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ {E}_p=M{\psi }_c^x+{\displaystyle \frac{Mx{\psi }_g^x}{a1}}\cdot ({\displaystyle \sum _{i=0}^N\frac{{(-1)}^i\cdot {({\psi }_s/{\psi }_g)}^{-i\lambda +x}}{-i\lambda +x}}-{\displaystyle \sum _{i=0}^N\frac{{(-1)}^i\cdot {({\psi }_c/{\psi }_g)}^{-i\lambda +x}}{-i\lambda +x}}) \end{gathered}$$

(16)

$$\begin{gathered} \mathrm{If} {\psi }_g>{\psi }_c\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ {E}_p=M{\psi }_c^x+{\displaystyle \frac{Mx{\psi }_g^x}{a1}}\cdot ({\displaystyle \sum _{i=0}^N\frac{{(-1)}^i\cdot {({\psi }_s/{\psi }_g)}^{(i+1)\lambda +x}}{(i+1)\lambda +x}}-{\displaystyle \sum _{i=0}^N\frac{{(-1)}^i\cdot {({\psi }_c/{\psi }_g)}^{(i+1)\lambda +x}}{(i+1)\lambda +x}}) \\ {\psi }_c<{\psi }_s<={\psi }_g\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \end{gathered}$$

(17)

$$\begin{gathered} E_p=E_p({\psi }_g)+{\displaystyle \frac{Mx{\psi }_g^x}{a1}}\cdot ({\displaystyle \sum _{i=0}^N\frac{{(-1)}^i\cdot {({\psi }_s/{\psi }_g)}^{-i\lambda +x}}{-i\lambda +x}}-{\displaystyle \sum _{i=0}^N\frac{{(-1)}^i}{-i\lambda +x}}) \\ {\psi }_s>{\psi }_g\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \end{gathered}$$

(18)

The value of $M$ is simulated randomly, and the value of $x$ is calculated based on the concept of diving evaporation. Let $z({\psi }_s)=D$ denote the depth of groundwater. The absolute value of the matric potential ${\psi }_s$, corresponding to $z({\psi }_s)$, can be determined using a trial algorithm. The variable ${\psi }_s$ is then substituted into the analytic solution for evaporation capacity, yielding the variable $E_p$. The optimal parameter group $(M,x)$ is selected by minimizing the error between the calculated evaporation capacity and the actual evaporation capacity. This parameter group is used for the power exponential function.

By employing $z({\psi }_s)$, $E_p$ and $E_g$ as constraints within the model, it is possible to determine the parameters $(M,x)$ of the power exponential function. Once these parameters are established, they can be utilized to calculate the vertical distribution of soil moisture, providing a detailed profile of moisture content throughout the soil layers.

2.3.2. The Analytic Solution of the Storage Capacity About Topographic Index Under Unsteady-State Conditions

According to Formulas (5)–(7), the storage capacity can be expressed as

$$\begin{gathered} W_{han}={\displaystyle {\int }_{z({\psi }_c)}^{z({\psi }_s)}({\theta }_r+({\theta }_s-{\theta }_r){(\frac{\psi }{{\psi }_b})}^{-\lambda })dz} (\mathrm{According} \mathrm{to} \mathrm{Formula} (5): {\psi }_m{\displaystyle \frac{1}{1+{(\psi /{\psi }_m)}^{p+x}}}d(\psi /{\psi }_m)=dz \\ ={\theta }_r(z({\psi }_s)-z({\psi }_c))+({\theta }_s-{\theta }_r){\displaystyle {\int }_{{\psi }_c}^{{\psi }_s}{(\frac{\psi }{{\psi }_b})}^{-\lambda }{\psi }_m(\frac{1}{1+{(\psi /{\psi }_m)}^{p+x}})d(\frac{\psi }{{\psi }_m})}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ ={\theta }_r(z({\psi }_s)-z({\psi }_c))+({\theta }_s-{\theta }_r){\psi }_b^{\lambda }{\psi }_m^{1-\lambda }{\displaystyle {\int }_{y_c}^{y_s}\frac{y^{-\lambda }}{1+y^{p+x}}dy\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}} \end{gathered}$$

(19)

$$\begin{gathered} W_f={\theta }_f\cdot (z({\psi }_s)-z({\psi }_c))-W_{han}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ =({\theta }_f-{\theta }_r)\cdot (z({\psi }_s)-z({\psi }_c))-({\theta }_s-{\theta }_r){\psi }_m^{1-\lambda }{\psi }_b^{\lambda }{\displaystyle {\int }_{y_c}^{y_s}\frac{y^{-\lambda }}{1+y^{p+x}}dy} \end{gathered}$$

(20)

$$\begin{gathered} \mathrm{when} y<1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ {\displaystyle \frac{y^{-\lambda }}{1+y^{p+x}}}={\displaystyle \sum _{i=0}^N{(-1)}^i}{y}^{i(p+x)-\lambda }, {\displaystyle \int \frac{y^{-\lambda }}{1+y^{p+x}}}dy={\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{y}^{i(p+x)+1-\lambda }}{i(p+x)+1-\lambda }} \end{gathered}$$

(21)

$$\begin{gathered} \mathrm{when} y>1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ {\displaystyle \frac{y^{-\lambda -p-x}}{y^{-p-x}+1}}={\displaystyle \sum _{i=0}^N{(-1)}^i}{y}^{-(i+1)(p+x)-\lambda }, {\displaystyle \int \frac{y^{-\lambda -p-x}}{y^{-p-x}+1}}dy={\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{y}^{-(i+1)(p+x)+1-\lambda }}{-(i+1)(p+x)+1-\lambda }} \end{gathered}$$

(22)

$$\begin{gathered} {\psi }_m<={\psi }_c, W_f=({\theta }_f-{\theta }_r)\cdot (z({\psi }_s)-z({\psi }_c))-({\theta }_s-{\theta }_r){\psi }_m^{1-\lambda }{\psi }_b^{\lambda }\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ ({\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{({\psi }_s/{\psi }_m)}^{-(i+1)(p+x)+1-\lambda }}{-(i+1)(p+x)+1-\lambda }-{\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{({\psi }_c/{\psi }_m)}^{-(i+1)(p+x)+1-\lambda }}{-(i+1)(p+x)+1-\lambda }}}) \end{gathered}$$

(23)

$$\begin{gathered} {\psi }_m>{\psi }_c, W_f=({\theta }_f-{\theta }_r)\cdot (z({\psi }_s)-z({\psi }_c))-({\theta }_s-{\theta }_r){\psi }_m^{1-\lambda }{\psi }_b^{\lambda }\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ ({\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{({\psi }_s/{\psi }_m)}^{i(p+x)+1-\lambda }}{i(p+x)+1-\lambda }-{\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{({\psi }_c/{\psi }_m)}^{i(p+x)+1-\lambda }}{i(p+x)+1-\lambda }}}) \\ {\psi }_c<{\psi }_s<={\psi }_m \end{gathered}$$

(24)

$$\begin{gathered} W_f=W_f({\psi }_m)+({\theta }_f-{\theta }_r)\cdot (z({\psi }_s)-z({\psi }_m))-({\theta }_s-{\theta }_r){\psi }_m^{1-\lambda }{\psi }_b^{\lambda }\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \\ ({\displaystyle \sum _{i=0}^N\frac{{(-1)}^i{({\psi }_s/{\psi }_m)}^{-(i+1)(p+x)+1-\lambda }}{-(i+1)(p+x)+1-\lambda }-{\displaystyle \sum _{i=0}^N\frac{{(-1)}^i}{-(i+1)(p+x)+1-\lambda }}}){\psi }_s>{\psi }_m \end{gathered}$$

(25)

$$z({\psi }_s)=z({\psi }_c+{\psi }_1)+sz\cdot {\displaystyle \frac{\max (t_{pi})-t_{pi}}{\max (t_{pi})-\min (t_{pi})}}$$

(26)

where $y={\displaystyle \frac{\psi }{{\psi }_m}}$; $y_c={\displaystyle \frac{{\psi }_c}{{\psi }_m}}$; $y_s={\displaystyle \frac{{\psi }_s}{{\psi }_m}}$; $z(\psi )$ signifies the vertical distance between the groundwater table and the point where the absolute value of the matric potential is $\psi$(mm); $z({\psi }_c)$ signifies the vertical distance between the groundwater table and the point where the absolute value of the matric potential is ${\psi }_c$ (mm); $W_{han}$ denotes the soil moisture content above the point where the absolute value of the matric potential is ${\psi }_c$ (mm); $W_f$ denotes the storage capacity(mm); $z({\psi }_c+{\psi }_1)$ represents the maximum groundwater depth of the wettest grid point, where the surface matric potential’s absolute value is ${\psi }_c+{\psi }_1$ and the storage capacity is $W_0$ (mm); $z({\psi }_c+{\psi }_1)+sz$ represents the maximum groundwater depth of the driest grid point (mm); $\max (t_{pi})$ denotes the maximum topographic index of the watershed; $\min (t_{pi})$ denotes the minimum topographic index of the watershed; and $t_{pi}$ denotes the topographic index of a specific grid point.

The evaporation capacity is assumed to be constant when determining the storage capacity (utilizing the evaporation capacity under the driest conditions). Subsequently, the variable $z({\psi }_s)$ is derived using Formula (26), which incorporates the topographic index of the grid point. The storage capacity is then determined through the application of Formulas (19)–(25). By mapping the distribution of topographic indices, a storage capacity curve under unsteady-state conditions is established.

2.3.3. Runoff Generation and Confluence Calculation

In this study, the distributed Xin’anjiang Model is employed to compute the dynamics of runoff generation and its subsequent confluence. By deriving the storage capacity curve under unsteady-state conditions, the model facilitates the acquisition of detailed insights into the temporal and spatial distribution of soil moisture, as well as the discharge patterns across the watershed.

2.3.4. The Simulation of the Storage Capacity

To validate the precision of the analytic solution for the storage capacity under unsteady-state conditions, four distinct simulation methods are employed:

• The analytic solution of the storage capacity under unsteady-state conditions.
• The numerical solution under unsteady-state conditions, where a very small time step is utilized to approximate the true value closely.
• The analytic solution of the storage capacity under equilibrium-state conditions, assuming all vertical evaporation fluxes are zero.
• The analytic solution of the storage capacity under steady-state conditions, where all vertical evaporation fluxes are set to the actual evaporation.

Unsteady-state conditions:

$$K(\psi )({\displaystyle \frac{\partial \psi }{\partial z}}-1)=M{\psi }^x$$

(27)

Equilibrium-state conditions:

$$K(\psi )({\displaystyle \frac{\partial \psi }{\partial z}}-1)=0$$

(28)

Steady-state conditions:

$$K(\psi )({\displaystyle \frac{\partial \psi }{\partial z}}-1)=Es$$

(29)

where $Es$ signifies the actual evaporation (mm/day) and $K(\psi )({\displaystyle \frac{\partial \psi }{\partial z}}-1)$ denotes the vertical evaporation flux.

By comparing the simulation results of the storage capacity from the numerical solution (which is considered nearly the true value) with the other three methods across various groundwater depths, the accuracy of these methods can be assessed. The precision of each method is determined by the degree to which its simulation results align with those of the numerical solution. A smaller difference between the storage capacity of the numerical solution and the result of another method indicates higher precision for that method.

2.3.5. Model Parameters and Evaluation Metrics

The distributed Xin’anjiang Model, which employs the analytic solution of the storage capacity under unsteady-state conditions, has been developed. The model integrates multiple parameters, including the difference between the maximum depth of the driest grid point and the maximum depth of the wettest grid point (SZ), soil parameters (l, also known as $\lambda$ and fyb, also known as ${\psi }_b$), difference between the saturated moisture content and the wilting moisture content (oo), evaporation conversion coefficient (eta), maximum storage capacity of free water storage reservoir at a specific location within the watershed (sm), outflow coefficient of interflow runoff (ki), outflow coefficient of groundwater runoff (kg), receding water coefficient of interflow runoff (ei), receding water coefficient of groundwater runoff (eg), confluence coefficient of the first day of surface runoff (es1), th parameter concerning soil evaporation (a1), saturated hydraulic conductivity coefficient (Ks), extreme depth of diving evaporation (dmax), upper tension water capacity (wum), lower tension water capacity (wlm), deep evaporation coefficient (ccc), parameter of diving evaporation (nn), and the difference between the absolute value of the surface matric potential of the wettest grid point when its depth is at its maximum and the value of ${\psi }_c$ (fy1, also known as ${\psi }_1$). All these parameters of the new model presented in this paper have been calibrated using the Shuffled Complex Evolution–Uncertainty Analysis (SCE-UA) algorithm [26].

Evaluating the simulation results and computational speed is essential and should be based on multiple objective functions. When assessing computational speed, this study considers two key metrics: the amount of computation required for a single stochastic parametric simulation and the time taken to complete 100 such simulations. Ideally, the larger the amount is and the shorter the time is, the faster the computational speed of the model will be. To assess the simulation accuracy of soil moisture, this study employs three evaluation indicators: the Nash efficiency coefficient (NSC), the root-mean-squared error (RMSE), and the determination coefficient (R²). For evaluating the simulation precision of discharge, two indices are used: NSC and RMSE. A higher NSC value indicates greater simulation precision for discharge or soil moisture. A higher R² value signifies greater simulation accuracy for soil moisture. A lower RMSE value suggests a better simulation effect for soil moisture or discharge. These indicators provide a comprehensive assessment of the model’s performance in simulating both soil moisture and discharge.

We conducted a comparative analysis of the objective functions across three models to assess their computational efficiency and accuracy. The Storage Capacity Model, despite its simplicity and reliance on readily available data, offers a modest calculation load and incorporates a basic physical mechanism. However, it oversimplifies the vertical hydrodynamic properties of the unsaturated zone and neglects the impact of setting off, leading to potential inaccuracies in soil moisture calculations. On the other hand, the DHSVM boasts a robust physical mechanism but demands extensive data and entails substantial computational demands. The accuracy of the DHSVM is compromised by data scarcity. The new model presented in this paper amalgamates the strengths of both predecessors while mitigating their respective weaknesses. It strikes a balance between calculation accuracy, computational speed, and the incorporation of a comprehensive physical mechanism. Our comparative analysis reveals that the new model outperforms the others in achieving this balance.

3. Results

3.1. The Comparison of the Storage Capacity Among Four Methods

To validate the accuracy of the analytic solution for the storage capacity under unsteady-state conditions, four simulation methods are utilized:

• Unsteady-State Solution: This is the analytic solution of the storage capacity under unsteady-state conditions.
• Numerical Solution for Unsteady State: The truncation error of this numerical solution is related to the time step. When the time step is set to 1 s, the truncation error is virtually zero, making the calculated value approximately the true value. The vertical distribution of soil moisture is calculated using this numerical solution, which takes into account the diving evaporation flux at the lower boundary and the actual evaporation flux at the upper boundary. The time step for these calculations is 1 s and the calculated value is almost the true value. To assess the accuracy of other methods, we compare their calculated values with those obtained from the numerical solution. A smaller difference between the values indicates higher accuracy for the method being compared.
• Equilibrium-State Solution: This is the analytic solution of the storage capacity under equilibrium-state conditions, where all vertical evaporation fluxes are zero.
• Steady-State Solution: This is the analytic solution of the storage capacity under steady-state conditions, where all vertical evaporation fluxes are set to the actual evaporation.

For this study, a typical soil is selected with the following parameters: ${\theta }_s-{\theta }_r=0.32,{\psi }_b=$$0.3 \mathrm{m},\lambda =0.4,p=3.2,Ks=1/{10}^6 \mathrm{m}/\mathrm{s}, a1=0.2$. The parameters for diving evaporation are set as $E_p=2 \mathrm{m}\mathrm{m}/\mathrm{d}\mathrm{a}\mathrm{y},d_{\max }=3.5 \mathrm{m},nn=2$. (In the new model presented in this paper, several soil-related parameters are considered: the saturated moisture content ${\theta }_s$, the wilting moisture content ${\theta }_r$, the intake potential ${\psi }_b$, the parameters specific to the Brooks–Corey model about soil pore size and hydraulic conductivity $\lambda ,p$, the saturated hydraulic conductivity coefficient $Ks$, and the soil evaporation parameter $a1$. Additionally, parameters related to diving evaporation are included, such as the extreme depth of deep evaporation $d_{\max }$ and the diving evaporation parameter $nn$). The vertical distributions of soil moisture are simulated for groundwater depths of 0.8 m, 1.3 m, and 1.8 m using the four methods described above. The accuracy of other methods is compared by the difference between the calculated values of the numerical solution and those of another method. The relationships between evaporative flux and matric potential for groundwater depths of 0.8 m, 1.3 m, and 1.8 m are analyzed for the unsteady-state solution and the numerical solution for unsteady state.

The vertical distributions of soil moisture and the relationships between the evaporative flux and matric potential for groundwater depths of 0.8 m, 1.3 m, and 1.8 m are as follows.

Figure 3, Figure 4 and Figure 5 reveal that the simulated vertical distributions of soil moisture produced by the unsteady-state solution and the numerical solution under unsteady-state conditions are nearly identical. When the groundwater depth is small, the simulated storage capacities from all four methods show minimal differences. However, as the depth increases, the discrepancies among the numerical solution, the equilibrium-state solution, and the steady-state solution become more pronounced. Specifically, at a depth of 1.8 m, the simulated storage capacities from the unsteady-state solution, the numerical solution, the equilibrium-state solution, and the steady-state solution are $131.11 \mathrm{mm},130.36 \mathrm{mm},101.86 \mathrm{mm}$ and $160.03 \mathrm{mm}$, respectively. The relative errors for the equilibrium-state solution and the steady-state solution are $-21.86\%$ and $22.76\%$, respectively. Compared to the $0.56\%$ relative error of the unsteady-state solution, these errors are significantly larger. Additionally, the actual evaporation calculated by the unsteady-state solution and the numerical solution at a depth of 1.8 m are $1.339$ $\mathrm{m}\mathrm{m}/\mathrm{d}\mathrm{a}\mathrm{y}$ and $1.344 \mathrm{mm}/\mathrm{day}$, respectively, with a relative error of only $-0.40\%$. The calculated diving evaporation at this depth is $0.47 \mathrm{mm}/\mathrm{day}$. (In the context of the numerical solution for the unsteady state, the truncation error is inherently linked to the time step size. When the time step is set to 1 s, the truncation error is negligible, effectively rendering the calculated actual evaporation equivalent to the true value. Subsequently, we employ the unsteady-state solution to estimate the actual evaporation. By comparing the calculated actual evaporation from the unsteady-state solution with the numerical solution for the unsteady state, we can ascertain the accuracy of the actual evaporation as determined by the unsteady-state solution).

At a groundwater depth of 1.3 m, the storage capacities calculated by the unsteady-state solution and the numerical solution are $59.45 \mathrm{mm}$ and $59.43 \mathrm{mm}$, respectively. The relative error between these two values is an insignificant $0.03\%$, indicating a high degree of consistency between the methods. Furthermore, the actual evaporation calculated by the unsteady-state solution and the numerical solution are $1.669 \mathrm{mm}/\mathrm{day}$ and $1.674 \mathrm{mm}/\mathrm{day}$, respectively. The relative error here is a minimal $-0.31\%$, suggesting that both methods provide a very similar and accurate estimation of actual evaporation. Additionally, the diving evaporation calculated at this depth is $0.80 \mathrm{mm}/\mathrm{day}$.

At a groundwater depth of 0.8 m, the storage capacities calculated by the unsteady-state solution and the numerical solution are $12.08 \mathrm{mm}$ and $12.20 \mathrm{mm}$, respectively. The relative error between these two values is a minor $-1.01\%$, indicating a high level of agreement between the methods. In addition, the actual evaporation calculated by the unsteady-state solution and the numerical solution are $1.938 \mathrm{mm}/\mathrm{day}$ and $1.914 \mathrm{mm}/\mathrm{day}$, respectively. The relative error here is a modest $1.25\%$, suggesting that both methods provide a reasonably close estimation of evaporation. Furthermore, the diving evaporation calculated at this depth is $1.20 \mathrm{mm}/\mathrm{day}$.

Figure 6, Figure 7 and Figure 8 reveal that the determination coefficient (R²) of the power exponential function trend line, which models the relationship between the evaporation flux and the absolute value of the matric potential using the numerical solution, consistently exceeds 0.94. This high R² value indicates that the power exponential function is a reliable choice for representing this relationship. Moreover, a comparison of the power exponential curves calibrated from the unsteady-state solution in this study with the trend lines of the numerical solution shows an extremely close similarity. (For the power exponential curves, when $D=0.8 \mathrm{m},M=2.03,x=$$0.437$. When $D=1.3 \mathrm{m}$, $M=1.135,x=$$0.291$. When $D=1.8 \mathrm{m}$, $M=0.631,x=$$0.313$. For the trend lines, when $D=0.8 \mathrm{m},M=$$1.902,x=$$0.447$. When $D=1.3 \mathrm{m}$, $M=1.116,x=$$0.392$. When $D=$$1.8 \mathrm{m}$, $M=$$0.605,x=$$0.339$). The phreatic evaporation and the actual evaporation derived from both the unsteady-state solution and the numerical solution for groundwater depths of 0.8 m, 1.3 m, and 1.8 m are found to be nearly identical. This observation is further supported by the visual analysis presented in Figure 6, Figure 7 and Figure 8.

3.2. Model Performance Evaluation

For the purpose of calibrating the model, meteorological and hydrological observation data from Tarrawarra Watershed, spanning from 12 December 1995, to 17 November 1997, have been selected. (Due to limitations in data availability, our study is based on meteorological, hydrological, and soil moisture data collected over a 2-year period. Given that our primary objective is to validate the accuracy of soil moisture measurements at 17 grid points, this duration of data is deemed sufficient for our analysis.) Throughout this observation period, the average daily precipitation recorded was 1.72 mm, with a maximum daily precipitation of 40 mm. The average daily discharge was 1.01 mm, peaking at a maximum daily discharge of 12.233 mm. The potential evapotranspiration was determined using the Penman–Monteith Equation [27,28], which was applied to the meteorological observation data. Following this, the parameters for the distributed Xin’anjiang Model, which employs the analytic solution of the storage capacity under unsteady-state conditions (referred to as the new model presented in this paper), the DHSVM, and the Storage Capacity Model were calibrated. Subsequently, the discharge and soil moisture for all three models were simulated. To assess the simulation effectiveness of the new model concerning discharge and soil moisture, a comparative analysis was conducted. This analysis involved evaluating the simulated discharge and soil moisture deficit within the 0–60 cm depth below the ground surface at 17 grid points from three different models against the corresponding observed data.

The soil moisture data from Tarrawarra Watershed pertains to the moisture levels at specific depths—15 cm, 30 cm, 45 cm, and 60 cm—below the ground surface across 17 grid points within the watershed. This information can be utilized to calculate the soil moisture deficit for each grid point

$$W_d=({\theta }_f-{\theta }_{15})\cdot 150+({\theta }_f-{\theta }_{30})\cdot 150+({\theta }_f-{\theta }_{45})\cdot 150+({\theta }_f-{\theta }_{60})\cdot 150$$

(30)

where $W_d$ signifies the soil moisture deficit of the grid point (mm); ${\theta }_{15},{\theta }_{30},{\theta }_{45},{\theta }_{60}$ denote the soil moisture at depths of 15 cm, 30 cm, 45 cm and 60 cm below the ground surface; and ${\theta }_f$ represents the field moisture capacity, in this paper, ${\theta }_f=0.38$.

In this study, the calibration of the parameters for the three models is performed using the SCE-UA algorithm. The calibration process is guided by an objective function that is defined as the sum of the Nash efficiency coefficients (NSCs) for both discharge and soil moisture deficit.

A comparison of the simulated discharge from three models against observed discharge data reveals the following. The distributed Xin’anjiang Model, which utilizes the analytic solution of the storage capacity under unsteady-state conditions (referred to as the new model presented in this paper), demonstrates the best simulation performance. It achieves a Nash efficiency coefficient (NSC) of 0.83 and a root-mean-squared error (RMSE) of 0.69 mm. The Storage Capacity Model exhibits a moderate simulation performance, with an NSC of 0.77 and an RMSE of 0.78 mm. The DHSVM, however, performs the least effectively among the three, with an NSC of 0.74 and an RMSE of 0.86 mm. The simulated discharge processes of all three models, alongside the observed discharge process, are depicted in Figure 9.

In this study, we employed the Sequential Chaotic Evolution with Univariate Marginal Distribution (SCE-UA) algorithm to calibrate model parameters. Initially, within the specified range of parameters, we randomly generated S sample points. These points served as the starting population for the algorithm. The SCE-UA algorithm then iteratively evolved these sample points (according to the rules of SCE-UA algorithm). At each iteration, new sample points were generated, and the objective functions for these new points were evaluated. The algorithm used these new points to replace the least fit individuals in the original sample set, thereby optimizing the objective functions. The evolutionary process continued until the convergence criteria were met. Specifically, the calculation halted when the improvement in the best objective function value fell below a predefined threshold or after a maximum number of iterations was reached. Upon convergence, the sample point with the optimal objective function value was identified. The parameters associated with this point were selected as the calibrated parameters for the model.

The parameters of the new model calibrated using the SCE-UA algorithm are presented in

Table 1. The calibrated parameters of the new model.

Parameter	sz	l	fyb	oo	eta	sm	ki
Value	1947.37	0.41	306.22	0.32	0.81	24.73	0.47
Range	1500–2500	0.1–0.8	150–500	0.2–0.4	0.7–1.2	10–50	0.1–0.7
Parameter	kg	ei	eg		es1	a1	Ks
Value	0.51	0.22	0.997		0.93	0.21	61.22
Range	0.1–0.7	0.5–0.95	0.99–0.999		0.5–1	0.1–0.3	10–100
Parameter	dmax	wum	wlm		ccc	nn	fy1
Value	3.52	21.53	53.91		0.23	2.03	324.57
Range	1–5	5–35	20–100		0.1–0.5	1–5	100–500

.

In the Tarrawarra Watershed, where each grid point measures 5 m in length and width with a total of 76 rows and 146 columns, the computational amount for a single stochastic parametric simulation of each of the three models was estimated. Concurrently, the time required for 100 stochastic parametric simulations for each model was recorded. The results of these evaluations are shown in

Table 2. The comparison of calculation speed among three models.

	The New Model	DHSVM	Storage Capacity Model
The calculation amount for a single stochastic parametric simulation	1,240,000	61,959,700	123,900
The calculation time required for 100 stochastic parametric simulations	2 min	100 min	12 s

.

Upon comparison, it was observed that the calculation speed of the new model presented in this paper is slightly slower than the Storage Capacity Model but significantly faster than the DHSVM.

The water deficits at 17 soil moisture observation points, as calculated using the parameters from the new model listed in Table 1, were analyzed across various time periods. The reliability of the simulation effect of the distributed Xin’anjiang Model, which employs the analytic solution of the storage capacity under unsteady-state conditions (as introduced in this paper), was assessed by comparing the simulated water deficits with the observed ones. The simulation results indicate that the Nash efficiency coefficient (NSC) for the water deficits at these 17 observation points ranged from 0.71 to 0.95. The root-mean-squared error (RMSE) was found to be below 21.11 mm. Moreover, the determination coefficient (R²) for 88% of the soil moisture observation points exceeded 0.90.

The changing process of simulated and observed soil moisture deficit within 0–60 cm below the ground surface at point S2 is depicted in Figure 10. The NSC, RMSE, and R² for the 17 simulated data points, as derived from three distinct models, are presented in

Table 3. The simulation results of water deficits within 0–60 cm below the ground surface at 17 soil moisture observation points for three models.

Point	NSC			RMSE			R2
	The New Model	Storage Capacity Model	DHSVM	The New Model	Storage Capacity Model	DHSVM	The New Model	Storage Capacity Model	DHSVM
S1	0.93	0.92	0.92	9.37	9.88	9.87	0.95	0.94	0.91
S2	0.92	0.66	0.68	8.43	17.03	16.47	0.93	0.74	0.84
S3	0.86	0.83	0.88	12.04	12.94	10.88	0.94	0.92	0.91
S4	0.84	0.83	0.90	12.90	13.05	10.11	0.92	0.91	0.93
S5	0.95	0.94	0.88	8.68	9.63	13.46	0.97	0.96	0.93
S6	0.93	0.94	0.85	10.04	8.67	14.29	0.92	0.96	0.82
S7	0.92	0.90	0.92	9.02	9.60	8.44	0.94	0.93	0.94
S8	0.92	0.90	0.86	10.62	11.79	14.42	0.94	0.92	0.89
S9	0.71	0.55	0.66	15.12	18.99	16.48	0.89	0.86	0.86
S10	0.82	0.81	0.79	18.61	19.18	20.14	0.96	0.94	0.92
S11	0.93	0.91	0.91	10.39	11.37	11.12	0.92	0.90	0.93
S12	0.90	0.88	0.84	10.62	11.57	13.62	0.94	0.93	0.91
S13	0.72	0.59	0.71	15.05	18.25	15.37	0.89	0.87	0.88
S14	0.93	0.92	0.91	10.23	10.90	11.81	0.92	0.91	0.91
S15	0.93	0.75	0.86	8.60	16.54	12.52	0.94	0.80	0.87
S16	0.79	0.76	0.75	21.11	21.89	22.46	0.93	0.91	0.90
S17	0.92	0.89	0.85	10.80	11.56	14.03	0.91	0.89	0.85
Average	0.88	0.82	0.83	11.86	13.70	13.83	0.93	0.90	0.89

.

A comparison among the simulation results of water deficits at 17 soil moisture observation points and discharge for the new model presented in this paper, the DHSVM, and the Storage Capacity Model reveals that the distributed Xin’anjiang Model, which utilizes the analytic solution of the storage capacity under unsteady-state conditions, demonstrates the best performance. Specifically, the Nash efficiency coefficient (NSC) for discharge is 0.83, and the root-mean-squared error (RMSE) is 0.69 mm. For water deficit, the Nash efficiency coefficient (NSC) ranges from 0.71 to 0.95, and the root-mean-squared error (RMSE) is below 21.11 mm. Moreover, the determination coefficient (R²) for 88% of the soil moisture observation points exceeds 0.90. The average absolute differences between the NSC of individual grid points and their respective average NSC values for the new model presented in this paper, the Storage Capacity Model, and the DHSVM are 0.0611, 0.0965, and 0.0694, respectively. Similarly, the average absolute differences between the RMSE of individual grid points and their respective average RMSE values are 2.78 mm, 3.51 mm, and 3.09 mm, respectively. Lastly, the average absolute differences between the R² of individual grid points and their respective average R² values are 0.0165, 0.0394, and 0.0311, respectively. These results indicate that the new model presented in this paper exhibits significantly smaller differences in NSC, RMSE, and R² compared to the other two models. Consequently, the new model demonstrates greater stability and more reliable performance across the grid points. Notably, for the wet points (S2, S9, S15), the new model’s simulation of water deficit is significantly improved due to the consideration of setting off effects. In contrast, the Storage Capacity Model assumes a zero-storage capacity for the wettest point, resulting in a notably smaller simulated water deficit. By adjusting the storage capacity of the wettest point to W₀ in this paper, the simulation effect for water deficits at the wet points has been markedly enhanced. In terms of computational speed, the new model presented in this paper is slightly slower than the Storage Capacity Model but substantially faster than the DHSVM.

4. Discussion

4.1. Advantages of the Model

In this study, we have developed the distributed Xin’anjiang Model, which incorporates the analytic solution of the storage capacity under unsteady-state conditions. This new model has been rigorously verified in the Tarrawarra Watershed of Australia, demonstrating the following advantages.

It is posited that the relationship between the evaporation flux and the absolute value of the matric potential follows a power exponential function. The parameters of this power exponential function are determined based on the groundwater depth, the evaporation capacity, and the phreatic evaporation. Subsequently, an analytical solution for the storage capacity under unsteady-state conditions is derived based on the established relationship. The solution is then merged with the distributed Xin’anjiang Model, facilitating the model’s distributed calculation process. This approach addresses the challenge of simplifying evaporation processes and the hydrodynamic characteristics of the unsaturated zone, thereby refining the model’s physical mechanism.

A comparative analysis was conducted among the unsteady-state solution, the numerical solution, the equilibrium-state solution, and the steady-state solution to evaluate their accuracy in predicting the storage capacity and the actual evaporation. The findings reveal that the relative errors for the storage capacity and the actual evaporation in the unsteady-state solution are minimal, not exceeding 1.01% and 1.25%, respectively. In contrast, the equilibrium- and steady-state solutions exhibit significantly higher relative errors of storage capacity, potentially exceeding 20% when the water depth is 1.8 m. This comparison underscores the superior precision of the unsteady-state solution, suggesting its reliability for modeling purposes. The unsteady-state solution effectively simulates both the storage capacity and the actual evaporation. The use of the power exponential function to model the relationship between the evaporation flux and the absolute value of matric potential is justified by two compelling arguments. Firstly, the coefficient of determination (R²) for the power exponential trend line of the numerical solution, correlating the evaporation flux with the absolute value of matric potential, exceeds 0.94, indicating a strong fit. Secondly, the power exponential curve calibrated based on the unsteady-state solution in this study closely mirrors the trend line of the numerical solution. Utilizing this relationship, the analytic solution presented here accurately predicts the vertical profile of soil moisture, aligning closely with observed data. This concordance further validates the efficacy of the unsteady-state solution in capturing the dynamics of vertical soil moisture distribution.

The simulation outcomes from three distinct models within the Tarrawarra Watershed provide clear evidence of the superiority of the new model presented in this paper. This model, which employs an analytical solution for the storage capacity under unsteady-state conditions, achieves a Nash efficiency coefficient (NSC) and a root-mean-square error (RMSE) for discharge of 0.83 and 0.69 mm, respectively. These results are notably superior to those of the Storage Capacity Model, which recorded 0.77 and 0.78 mm, and the DHSVM, with 0.74 and 0.86 mm. Furthermore, the new model’s performance in predicting soil moisture is also commendable, with an average NSC, RMSE, and R² of 0.88, 11.86 mm, and 0.93. These figures outperform the Storage Capacity Model’s 0.82, 13.70 mm, and 0.90, as well as the DHSVM’s 0.83, 13.83 mm, and 0.89. This indicates a more precise representation of soil moisture dynamics. In terms of computational efficiency, the new model demonstrates a competitive edge. It requires only 2 min to complete 100 stochastic parametric simulations, which, while slightly slower than the Storage Capacity Model’s 12 s, is significantly faster than the DHSVM’s 100 min. This suggests that the new model offers a balance between accuracy and computational speed. Additionally, the new model’s consideration of the impact of setting off has notably enhanced the simulation accuracy at wet points (S2, S9, S15). For instance, the NSC for S2 has improved from 0.66 (Storage Capacity Model) to 0.92 (new model), for S9 from 0.55 to 0.71, and for S15 from 0.75 to 0.93. This improvement underscores the model’s enhanced ability to capture the nuances of hydrological processes at specific locations. In summary, the new model presented in this paper effectively simulates the spatial and temporal variations of hydrological variables and the watershed’s discharge process. It also exhibits a commendable calculation speed and demonstrates an accurate simulation of water deficits at wet points, making it a valuable tool for hydrological modeling and analysis.

The data required for the new model, such as terrain and soil characteristics, are easily obtainable. Relative to preceding distributed Xin’anjiang Models, the novel model presented in this paper has been carefully crafted to achieve a balance between computational velocity, predictive accuracy, and compliance with physical principles [14,15,16,17]. Thanks to these enhancements, the applicability of the distributed Xin’anjiang Model has been significantly broadened. It now incorporates an analytical solution for the storage capacity under unsteady-state conditions, which further refines its performance. The new model is poised for application in various fields, including hydrology and water resource simulation, ecological impact assessment, and hydrological forecasting in regions with limited data. It is designed to address the limitations of existing distributed hydrological models, particularly those related to data scarcity and the substantial computational demands [9,11,12].

4.2. Limitations and Outlook

In the current model, the evaporation capacity is assumed to be constant when determining the storage capacity (utilizing the evaporation capacity under the driest conditions). However, it is acknowledged that the evaporation capacity is subject to continuous variation, which in turn affects the vertical distribution of soil moisture [10]. Therefore, future studies should incorporate this dynamic factor into the model to enhance its accuracy and applicability. Additionally, the model should account for the diversity of soil types and the impact of vegetation on the hydrological cycle [1,2,10,29]. These factors are crucial for a comprehensive understanding of watershed behavior. Future research and discussions should focus on integrating soil and vegetation dynamics into the model to provide a more nuanced and realistic representation of watershed processes. Furthermore, the current study is constrained by the use of limited data and the validation of the new model in a single watershed. Future research endeavors should prioritize the integration of additional data sets and the expansion of model validation to encompass a broader range of watersheds. This approach will not only enhance the robustness of the model but also facilitate a more comprehensive understanding of its applicability across different environmental conditions.

5. Conclusions

In this research, an analytical solution for the storage capacity under unsteady-state conditions has been derived and subsequently integrated into the distributed Xin’anjiang Model. This approach refines the model’s physical mechanism. A comparison is made between the unsteady-state solution, the numerical solution, the equilibrium-state solution, and the steady-state solution. It is observed that the power exponential curve calibrated from the unsteady-state solution closely resembles the curve simulated by the numerical solution. The unsteady-state solution presented in this study effectively simulates the vertical distribution of soil moisture and the actual evaporation. Furthermore, a comparative analysis of the simulation results from three models applied to the Tarrawarra Watershed reveals that the new model proposed in this paper provides the most accurate simulation of the spatial and temporal variations of soil moisture and the discharge process. Although the calculation speed of the new model is slightly slower than that of the Storage Capacity Model, it is significantly faster than the DHSVM. The improved simulation of wet points in the new model is attributed to the consideration of setting off effects, which results in a noticeably better performance compared to the Storage Capacity Model.

This study demonstrates that the newly introduced model achieves a commendable balance among computational efficiency, predictive accuracy, and adherence to physical principles. The model necessitates data inputs such as terrain and soil characteristics, which are typically readily available. This novel model holds promise for application across various sectors, offering a significant advancement over current distributed hydrological models. The model is capable of overcoming the limitations associated with the limited data availability and the high computational demands.

However, the current version of the new model lacks a specific method for considering the temporal variability in evaporation capacity. It also does not differentiate between soil types or account for the effects of vegetation. Moreover, the new model currently relies on a limited data set and has been validated exclusively within a single watershed. To enhance the distributed Xin’anjiang Model, future research should focus on bridging the identified gaps. This can be achieved by developing a method that accounts for the dynamic changes in evaporation capacity during soil moisture calculations. Additionally, the method should include the classification of various soil types, the integration of vegetation impacts within the model’s framework, and the incorporation of additional data sets. Furthermore, it is essential to validate the refined model across a broader range of watersheds to ensure its applicability and robustness.