Two-dimensional differential form of distributed Xinanjiang model

The distributed hydrologic models (DHMs) evolved from lumped hydrologic models, inheriting their modeling philosophy along with persistent numerical-error issues. Historically, these models tend to use established one-dimensional (1D) methods for slope concentration, which often struggle to effectively represent complex terrains. In this study, we formulated a purely differential form of mathematical equations for the distributed Xinanjiang model and developed a fully coupled numerical solution framework. We also introduced two-dimensional (2D) diffusion wave equations for surface slope concentration and derived 2D linear reservoir equations for subsurface slope concentration to replace their 1D counterparts. This culminated in the development of a two-dimensional differential form of the distributed Xinanjiang (TDD-XAJ) model. Two numerical experiments and the application of the TDD-XAJ model in a humid watershed were conducted to demonstrate the model. Our results suggested that (a) numerical errors in the existing distributed Xinanjiang model are significant and may be exacerbated by a potential terrain amplification effect, which could be effectively controlled by the fully coupled numerical framework within the TDD-XAJ model; (b) the 2D slope concentration methods showed enhanced terrain capture ability and eliminated the reliance on the flow direction algorithms used in 1D methods; and (c) the TDD-XAJ model exhibited improved simulation capabilities compared to the existing model when applied in the Tunxi watershed, particularly for flood volume. This study emphasizes the need to revisit DHMs which stem from lumped hydrological models, focusing on model equations and numerical implementations, which could enhance model performance and benefit the hydrological modeling community.

1 Introduction

The distributed hydrological model (DHM) is a sophisticated perceptual, mathematical, and computational framework that simulates the spatial and temporal dynamics of hydrological processes across watersheds, particularly addressing the rainfall–runoff partitioning (Fatichi et al., 2016; Kampf and Burges, 2007; Loritz et al., 2021; Paniconi and Putti, 2015). By effectively accounting for climate conditions and critical-zone characteristics, DHMs have been widely used for hypothesis testing (Clark et al., 2016; Nippgen et al., 2015), flood forecasting (Perrini et al., 2024), future scenario projections (Ul Hassan et al., 2024), and other applications. The development of DHMs benefits from advancements in observation methods and computational techniques, as well as crucial insights gained from hydrological experiments (Bisht and Riley, 2019; Dongarra and Keyes, 2024; Paniconi and Putti, 2015). The equations used in DHMs can be classified into two categories: rigorous mathematical–physical equations simplified from the Navier–Stokes equations and parameterized equations from the lumped hydrological models. The emphasis on mathematical–physical equations underscores their theoretical grounding, while parameterized equations highlight their consistency with observations (Beven, 2002; Freeze and Harlan, 1969; Reggiani et al., 1998). This integration of both approaches expands the available equation space for DHMs by absorbing the merits of lumped hydrological models (Beven, 2002; Tran et al., 2018).

The Xinanjiang (XAJ) model (Zhao, 1992) is widely regarded to be a standard hydrological model in China (Chen et al., 2023; Zhao et al., 2023) and has been recognized worldwide (He et al., 2023; Knoben et al., 2020; Taheri et al., 2023; Wang et al., 2023). The XAJ model assumes that runoff can only occur once the soil tension water capacity is satisfied, following a saturation–excess runoff generation mechanism (Nan et al., 2024; Zhao, 1992). A distinctive feature of the XAJ model is its use of the Pareto distribution to effectively represent the spatial distribution of soil tension water capacity within a watershed in a parsimonious manner (Taheri et al., 2023). The development of the XAJ model has evolved through three distinct phases over the last 60 years.

• Phase 1 (1963–1980). This phase marks the initial formulation of the XAJ model. It began in 1963 when the saturation–excess runoff mechanism was first introduced (Zhao and Zhuang, 1963) and concluded in 1980 when the model was formally named the Xinanjiang model (Zhao et al., 1980). Key contributions during this phase included the establishment of a runoff generation module based on probability distribution curves and a top-down procedure that separated the total runoff into distinct components – namely surface runoff and groundwater – before routing them to the channel.

• Phase 2 (1980–2002). The second phase focused on enhancing the structure and accuracy of the XAJ model. Three major improvements were made: (a) the evapotranspiration module was refined by dividing the soil horizon into three layers, addressing the issue of underestimating evapotranspiration after prolonged dry periods; (b) interflow was introduced as a new runoff source, inspired by the progress in hillslope hydrology during the International Hydrological Decade (McGuire et al., 2024); and (c) the original hydrograph method for slope and channel concentration was replaced with a linear reservoir and lag-routing method, leading to the formation of the widely used three-water-source lumped XAJ model (Zhao, 1992).

• Phase 3 (2002–present). The third phase is characterized by efforts to transition the XAJ model from a lumped to a distributed version, aligning it with other contemporary models like TOPMODEL and HBV (Beven et al., 2021; Seibert et al., 2022). It is noteworthy that the original lumped XAJ model has undergone continuous evolution alongside the distributed version of the model (Ouyang et al., 2025). Although initial efforts began earlier (Lu et al., 1996), 2002 was a notable milestone due to Beven's alternative blueprint (Beven, 2002). This blueprint emphasized the importance of observational consistency as a core requirement of the physically based property of DHMs, enabling the integration of parameterized equations from lumped hydrological models. The transformation involved the application of the runoff generation modules to smaller computational units (e.g., sub-basins or grids) and the implementation of distributed hydrological or hydraulic runoff concentration methods (Chen et al., 2024; Fang et al., 2017; Liu et al., 2009; Su et al., 2003). Grid-based models, which utilize structured grids for data input and output, are easier to preprocess, implement, and visualize (Shu et al., 2024). This approach fosters connections with other disciplines and benefits from ongoing community support (Chen et al., 2023; Wu et al., 2024; Zhang et al., 2024), making grid-based distributed models more popular and a focal point of this study. The Grid-XAJ (GXAJ) model is a representative achievement of the existing distributed Xinanjiang model (Yao et al., 2009; Yao et al., 2012).

A limitation of the current generation of the distributed XAJ model is its reliance on a one-dimensional slope concentration method. This method assumes that lateral-flow characteristics that deviate from the dominant flow direction can be neglected (Hong and Mostaghimi, 1997). Consequently, one-dimensional (1D) methods for channel concentration can be adapted to represent slope concentration, ranging from simpler approaches like linear reservoir or Muskingum methods to more complex options, such as kinematic and diffusion wave models (Clark et al., 2015; Paniconi and Putti, 2015; Todini, 2007). Typically, the 1D slope concentration method is combined with a single-flow-direction algorithm, such as the well-known eight-direction (D8) algorithm (O'Callaghan and Mark, 1984), to determine the dominant flow direction (Zhu et al., 2013). This practice is common in grid-based distributed hydrological models, such as WECOH (Nippgen et al., 2015) and HydroPy (Stacke and Hagemann, 2021). However, the 1D slope concentration method has long been criticized for its inability to accurately simulate water flow in complex terrains and its failure to adequately represent the effects of microtopography (Hong and Mostaghimi, 1997; Liu et al., 2004).

Another major deficiency of the current generation of the distributed XAJ model is the numerical-error issue inherited from its lumped counterparts, a common phenomenon in lumped models (Clark and Kavetski, 2010; Gupta et al., 2012). These lumped models typically rely on the difference form of equations formulated over discrete time periods using simple numerical techniques, such as the first-order explicit Euler method and operator splitting (Clark and Kavetski, 2010; Santos et al., 2018; Schoups et al., 2010; Woldegiorgis et al., 2023). Although more advanced methods, such as the MacCormack scheme (MacCormack, 1982), had been applied for runoff concentration, numerical errors persisted during the development of the distributed XAJ model. This ongoing issue can be attributed to the direct application of the difference form of runoff generation equations from the lumped XAJ model. Moreover, a loosely coupled numerical implementation framework had been developed, where the total number of different runoff sources was calculated separately using difference form equations, assuming their intensities remained constant over discrete time intervals when fed into the differential form of runoff concentration equations (Yao et al., 2012). This framework further complicated the existing numerical-error problem. Such numerical errors can deteriorate model performance (Clark and Kavetski, 2010; Woldegiorgis et al., 2023), complicate parameter calibration (Kavetski and Clark, 2010), and even impact the effectiveness of physically informed machine learning methods (Song et al., 2024). Furthermore, the effect of numerical errors tends to intensify with increased precipitation, leading to greater uncertainty in model applications (La Follette et al., 2021).

A differential form of the mathematical framework that describes the runoff generation and concentration process for the distributed XAJ model is required to address the numerical-error issue, which further enables a fully coupled numerical implementation framework. This effort is based on two key ideas. First, Zhao et al. (2023) introduced a differential form of the lumped XAJ model, which provides ordinary differential equations (ODEs) to describe the grid-scale runoff generation process. These ODEs can be integrated with partial differential equations (PDEs) and other ODEs that represent the runoff concentration process, resulting in complete mathematical equations of the distributed XAJ model. Second, we employed Qu's strategy for handling hydrological processes involving mixed PDEs and ODEs, utilizing the finite-volume method (FVM) to spatially discretize the PDEs (Qu and Duffy, 2007). This strategy resulted in a system of ODEs that can be solved in a unified manner.

The specific objectives of this study were as follows:

1. to introduce or derive 2D slope concentration methods for surface and subsurface runoff and to compare their performance with 1D methods;

2. to evaluate and control the numerical errors arising from the loosely coupled numerical implementation framework in the current generation of the distributed XAJ model;

3. to propose a new distributed XAJ model incorporating the differential form of equations, a fully coupled numerical implementation, and 2D slope concentration methods.

The remainder of this paper is organized as follows. Section 2 outlines the distributed XAJ model we developed. Section 3 introduces the numerical experiments and model application, including the details of the test cases and study area used. Section 4 discusses the results of the numerical experiments and model application. Section 5 provides the summary and conclusions of this study.

2 Model theory and development

2.1 General overview

The two-dimensional differential form of the distributed XAJ (TDD-XAJ) model inherits the merits of the original lumped XAJ model, incorporating 2 decades of efforts to transform it into a distributed model. The TDD-XAJ model uses a grid-based structure, assuming that meteorological forcing and model parameters are uniform within each grid but vary between grids. The key hydrological processes in the model are categorized into runoff generation, which is calculated at the grid level, and runoff concentration, which is calculated between grids. The soil column within each grid is divided into three layers for evapotranspiration calculation. Precipitation reaches the surface as net precipitation after accounting for actual evapotranspiration losses. Total runoff is calculated based on net precipitation using the saturation–excess runoff mechanism and is then partitioned into surface runoff, interflow, and groundwater runoff for further concentration. Water moves in two dimensions across the slope, with surface and subsurface runoff (interflow and groundwater) draining into the channels. The water in the channels is routed to the outlet in a one-dimensional flow pattern. Bidirectional water exchange between the slope surface and channel is considered under varying hydraulic conditions. Specifically, when the channel water surface elevation exceeds that of the slope surface, water flows from the channel back onto the slope surface. The diagram of the TDD-XAJ model is shown in Fig. 1.

Figure 1 Diagram of the two-dimensional differential form of the distributed Xinanjiang model.

The square digital elevation model (DEM) cells are used as fundamental computational units with a grid size of Δx (m). The channel is aligned either along the axis or diagonally through the center of the grid cell and is divided into segments corresponding to the grid cell with which it intersects. The river network consists of several channels, each assumed to have uniform hydraulic elements. Branching is not considered, meaning a channel may receive flow from multiple upstream channels and discharge into no more than one downstream channel. Only the length of the mainstream segment is considered, while tributary inflows are treated as being equivalent to the source term. The length of a channel segment Δl can be Δx, $(\mathrm{1}+\sqrt{\mathrm{2}})\mathrm{\Delta }x/\mathrm{2}$ or $\sqrt{\mathrm{2}}\mathrm{\Delta }x$ based on its flow direction. The distance between two neighboring grid centers δl is constrained to either Δx or $\sqrt{\mathrm{2}}\mathrm{\Delta }x$.

There are five key differences between the TDD-XAJ model and the GXAJ model, which is a representative of the existing distributed Xinanjiang model:

1. Revised module partitioning method. The evapotranspiration calculation module and the runoff generation module from the lumped XAJ model have been combined into a single module in the TDD-XAJ model as they share the same state variable. This adjustment benefits mathematical description and maintains consistency with previous work.

2. Differential-form runoff generation equations. Unlike existing models that use difference form equations derived with simple numerical techniques, the TDD-XAJ model adopts differential-form runoff generation equations for improved numerical accuracy.

3. Two-dimensional diffusion wave equations. Instead of the one-dimensional (1D) diffusion wave equations, the TDD-XAJ model employs two-dimensional (2D) diffusion wave equations to better represent surface water flow across the slope.

4. Extended subsurface runoff concentration equations. The 1D linear reservoir method previously used for subsurface runoff concentration has been extended to be 2D in the TDD-XAJ model, allowing for better representation of subsurface flow and aligning with the surface runoff concentration method for consistency.

5. Fully coupled numerical implementation framework. The TDD-XAJ model uses a fully coupled numerical approach, where the model's equations are assembled into a system of ODEs and are solved simultaneously. In contrast, the GXAJ model employed a loosely coupled method, with the difference form of the equations for runoff generation and the differential form of the equations for runoff concentration.

2.2 Mathematical structure

2.2.1 Evapotranspiration and runoff calculation

The actual evapotranspiration is calculated by applying the three-layer evapotranspiration (TLE) formulas of the lumped XAJ model. The TLE concept divides the soil horizon into upper, lower, and deep layers. Each layer has a tension water storage and storage capacity. Depending on the magnitude of precipitation intensity and potential evapotranspiration intensity, the tension water storages of these soil layers are replenished by precipitation or lost by evapotranspiration in a top-to-bottom sequence. When precipitation exceeds potential evapotranspiration intensity, the surplus (net precipitation) immediately enters the upper soil layer. If the tension water storage of the upper or lower soil layer exceeds its water storage capacity, the excess is transferred further to the lower or deep soil layer. Conversely, during conditions of less precipitation, evaporative losses are sequentially deducted from the soil layers, starting from the upper layer and moving to the deep layer. The governing differential equations of TLE formulas are written as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(\begin{array}{c}{W}_{\mathrm{u}}\\ W_{\mathrm{l}}\\ W_{\mathrm{d}}\end{array}\right)=\left(\begin{array}{c}{P}_{\mathrm{n}}-R-E_{\mathrm{u}}-I_{\mathrm{u}}\\ I_{\mathrm{u}}-E_{\mathrm{l}}-I_{\mathrm{l}}\\ I_{\mathrm{l}}-E_{\mathrm{d}}\end{array}\right),\end{array}$$

where t is time (s), and W_u, W_l, and W_d are three state variables which represent the tension water storage of the upper, lower, and deep soil layers (mm). The variables on the right-hand side of Eq. (1) are all model fluxes with units of mm s⁻¹, which include net precipitation intensity (P_n); total runoff intensity (R); actual evapotranspiration intensity from the upper (E_u), lower (E_l), and deep (E_d) soil layers; recharge intensity from the upper (I_u) and lower (I_l) soil layers to next soil layer when the tension water capacity is satisfied.

The tension water storage capacity curve (TWSCC) is used to calculate R, which depicts the influence of the spatial heterogeneity of tension water storage capacity on the runoff generation process using the Pareto distribution. The constitutive equations of R, which consider the effect of the impervious area, are written as follows:

$$\begin{array}{}\text{(2)}& {\displaystyle }\begin{array}{rl}{f}_{\mathrm{w}}& =A_{\text{imp}}+\\ & (\mathrm{1}-A_{\text{imp}})\left[\mathrm{1}-{\left(\mathrm{1}-{\displaystyle \frac{W_{\mathrm{u}}+W_{\mathrm{l}}+W_{\mathrm{d}}}{W_{\text{um}}+W_{\text{lm}}+W_{\text{dm}}}}\right)}^{b/(\mathrm{1}+b)}\right],\end{array}\text{(3)}& {\displaystyle }R=P_{\mathrm{n}}{f}_{\mathrm{w}},\end{array}$$

where f_w denotes the runoff coefficient or the ratio of areas where tension water capacity is satisfied (–); A_imp is the ratio of the impervious area (–); b is the TWSCC exponent (–); and W_um, W_lm, and W_dm are the tension water storage capacity of the three soil layers (mm).

The constitutive equations of the remainder of the model fluxes in Eq. (1) are given as follows:

$$\begin{array}{}\text{(4)}& \left(\begin{array}{c}{P}_{\mathrm{n}}\\ E_{\mathrm{n}}\\ E_{\mathrm{u}}\\ I_{\mathrm{u}}\\ E_{\mathrm{l}}\\ I_{\mathrm{l}}\\ E_{\mathrm{d}}\end{array}\right)=\left(\begin{array}{c}max(P_{\text{obs}}-K_{\mathrm{e}}{E}_{\text{obs}},\mathrm{0})\\ max(K_{\mathrm{e}}{E}_{\text{obs}}-P_{\text{obs}},\mathrm{0})\\ max(W_{\mathrm{u}},\mathrm{0})\cdot E_{\mathrm{n}}\\ max\left[\text{sgn}\right(W_{\mathrm{u}}-W_{\text{um}}),\mathrm{0}]\cdot (P_{\mathrm{n}}-R-E_{\mathrm{u}})\\ max(W_{\mathrm{l}},\mathrm{0})\cdot max(c,W_{\mathrm{l}}/W_{\text{lm}})\cdot (E_{\mathrm{n}}-E_{\mathrm{u}})\\ max\left[\text{sgn}\right(W_{\mathrm{l}}-W_{\text{lm}}),\mathrm{0}]\cdot (I_{\mathrm{u}}-E_{\mathrm{l}})\\ max(W_{\mathrm{d}},\mathrm{0})\cdot max\left[c\right(E_{\mathrm{n}}-E_{\mathrm{u}})-E_{\mathrm{l}},\mathrm{0}]\end{array}\right),\end{array}$$

where P_obs and E_obs are observed precipitation and pan evaporation intensity (mm s⁻¹), K_e is the coefficient of potential evapotranspiration to pan evaporation (–), E_n denotes the net evapotranspiration intensity (mm s⁻¹), and c is the coefficient of deep-soil-layer evapotranspiration (–).

2.2.2 Runoff separation

According to the top-down modeling philosophy of the XAJ model, the total runoff is calculated first and then separated into different runoff sources. A distinction is made between impervious and pervious areas when separating the total runoff; the total runoff on impervious areas is treated as surface runoff directly, while the total runoff on pervious areas is further divided into surface runoff, interflow, and groundwater runoff. The total amount of surface runoff is the summation of surface runoff from both the impervious and pervious areas. The separation procedure in pervious areas combines the free-water storage capacity curve (FWSCC) and the linear reservoir method. According to the Dunne saturation–excess runoff mechanism, it is assumed that the surface runoff only occurs when free-water storage capacity is satisfied. The spatial heterogeneity of free-water storage capacity is accounted for with the same Pareto distribution used in TWSCC. The interflow and groundwater runoff are assumed to show linear outflow depending on the free-water storage. The governing equation and constitutive equations of the runoff separation module are written as follows:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{S}_{\mathrm{0}}}{\mathrm{d}t}}=P_{\mathrm{n}}-(R_{\text{ps}}+R_{\mathrm{i}}+R_{\mathrm{g}})/(f_{\mathrm{w}}-A_{\text{imp}}),\text{(6)}& {\displaystyle }\begin{array}{rl}& \left(\begin{array}{c}{R}_{\text{ps}}\\ R_{\mathrm{i}}\\ R_{\mathrm{g}}\end{array}\right)=\\ & (f_{\mathrm{w}}-A_{\text{imp}})\left(\begin{array}{c}{P}_{\mathrm{n}}-P_{\mathrm{n}}(\mathrm{1}-S_{\mathrm{0}}/S_{\mathrm{m}}{)}^{\text{ex}/(\mathrm{1}+\text{ex})}\\ -K_{\mathrm{i}}{S}_{\mathrm{0}}\ln (\mathrm{1}-K_{\mathrm{i}}-K_{\mathrm{g}})/(K_{\mathrm{i}}\mathrm{\Delta }T+K_{\mathrm{g}}\mathrm{\Delta }T)\\ -K_{\mathrm{g}}{S}_{\mathrm{0}}\ln (\mathrm{1}-K_{\mathrm{i}}-K_{\mathrm{g}})/(K_{\mathrm{i}}\mathrm{\Delta }T+K_{\mathrm{g}}\mathrm{\Delta }T)\end{array}\right),\end{array}\text{(7)}& {\displaystyle }{R}_{\mathrm{s}}=P_{\mathrm{n}}{A}_{\text{imp}}+R_{\text{ps}},\end{array}$$

where S₀ is free-water storage (mm); S_m is free-water storage capacity (mm); R_ps, R_i, and R_g are surface runoff intensity from the pervious areas, interflow intensity, and groundwater intensity (mm s⁻¹); ex is the FWSCC exponent (–); K_i and K_g are the interflow and groundwater outflow coefficient (–); ΔT is the time interval of input forces (s); and R_s is total surface runoff intensity (mm s⁻¹).

2.2.3 Surface runoff concentration

The governing equations for the slope concentration process of the surface runoff are the two-dimensional diffusion wave equations (Gottardi and Venutelli, 2008), which are given as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial \mathit{U}}{\partial t}}+{\displaystyle \frac{\partial \mathit{F}\left(\mathit{U}\right)}{\partial x}}+{\displaystyle \frac{\partial \mathit{G}\left(\mathit{U}\right)}{\partial y}}=\mathit{S}\left(\mathit{U}\right),\text{(9)}& {\displaystyle }\begin{array}{rl}& \mathit{U}=\left[\begin{array}{c}{h}_{\mathrm{s}}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right]\mathit{F}=\left[\begin{array}{c}{h}_{\mathrm{s}}u\\ g{h}_{\mathrm{s}}^{\mathrm{2}}/\mathrm{2}\\ \mathrm{0}\end{array}\right]\mathit{G}=\left[\begin{array}{c}{h}_{\mathrm{s}}v\\ \mathrm{0}\\ g{h}_{\mathrm{s}}^{\mathrm{2}}/\mathrm{2}\end{array}\right]\\ & \mathit{S}=\left[\begin{array}{c}{\mathit{\varphi }}_{\mathrm{s}}\\ g{h}_{\mathrm{s}}(S_{\text{ox}}-S_{\text{fx}})\\ g{h}_{\mathrm{s}}(S_{\text{oy}}-S_{\text{fy}})\end{array}\right],\end{array}\end{array}$$

where h_s is the surface water depth (m); u and v are the surface flow velocity in the x and y directions, respectively (m s⁻¹); g denotes the acceleration due to gravity (m s⁻²) and is taken to be 9.81; S_ox and S_oy are the surface bottom slope term in the x and y directions, respectively (–); $S_{\text{ox}}=-\partial z_{\mathrm{s}}/\partial x$, and $S_{\text{oy}}=-\partial z_{\mathrm{s}}/\partial y$; z_s represents the surface elevation (m); S_fx and S_fy denote the surface friction term in the x and y directions, respectively (–); and ϕ_s is the source term (m s⁻¹).

The cell-centered finite-volume method (FVM) is used to spatially discretize the two-dimensional (2D) diffusion wave equations based on a square-structured mesh (Jain and Singh, 2005). The mesh layout and stencils used are shown in Fig. 2a.

Figure 2 Diagram of spatial discretization for runoff concentration equations.

Within the framework of the FVM, we integrate Eq. (8) over the control volume (Ω, red square in Fig. 2a) and introduce Green's theorem to transform the double integral over Ω into a line integral around its boundary Γ, which turns out to be

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial t}}{\int }_{\mathrm{\Omega }}\mathit{U}\mathrm{d}\mathrm{\Omega }+{\oint }_{\mathrm{\Gamma }}\mathit{F}\left(\mathit{U}\right)\cdot n\mathrm{d}l+{\oint }_{\mathrm{\Gamma }}\mathit{G}\left(\mathit{U}\right)& \cdot n\mathrm{d}l=\\ & {\int }_{\mathrm{\Omega }}\mathit{S}\left(\mathit{U}\right)\mathrm{d}\mathrm{\Omega }.\end{array}\end{array}$$

This is done with the incorporation of Eq. (9) and by further assuming that the surface water is evenly distributed in Ω. Equation (8) could be written in the ODE form after the spatial discretization of the FVM, which could be expressed as follows:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{h}_{\mathrm{s}}^{i,j}}{\mathrm{d}t}}=-& \left({\displaystyle \frac{h_{\mathrm{s},\text{eff}}^{i-\mathrm{1}/\mathrm{2},j}{u}^{i-\mathrm{1}/\mathrm{2},j}-h_{\mathrm{s},\text{eff}}^{i+\mathrm{1}/\mathrm{2},j}{u}^{i+\mathrm{1}/\mathrm{2},j}}{\mathrm{\Delta }x}}\right.\\ & \left.+{\displaystyle \frac{h_{\mathrm{s},\text{eff}}^{i,j-\mathrm{1}/\mathrm{2}}{v}^{i,j-\mathrm{1}/\mathrm{2}}-h_{\mathrm{s},\text{eff}}^{i,j+\mathrm{1}/\mathrm{2}}{v}^{i,j+\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta }x}}\right)\\ & +{\mathit{\varphi }}_{\mathrm{s}}^{i,j},\end{array}\end{array}$$

where i and j denote the index of Ω, $h_{\mathrm{s},\text{eff}}^{i+\mathrm{1}/\mathrm{2},j}$ and $h_{\mathrm{s},\text{eff}}^{i,j+\mathrm{1}/\mathrm{2}}$ are the effective water depth at the eastern and northern boundaries of Ω (m), $u^{i+\mathrm{1}/\mathrm{2},j}$ and $v^{i,j+\mathrm{1}/\mathrm{2}}$ are the surface flow velocity in the direction of the external normal vector at the eastern and northern boundaries of Ω (m s⁻¹), and Δx is the grid size (m). ϕ_s includes total surface runoff and the interaction discharge between the surface and channel Q_sc (m³ s⁻¹), whereby ${\mathit{\varphi }}_{\mathrm{s}}=R_{\mathrm{s}}/\mathrm{1000}-Q_{\text{sc}}/\mathrm{\Delta }{x}^{\mathrm{2}}$. The effective water depth at the boundary is determined according to the water depth and surface elevation on the left and right sides of the boundary (Bates et al., 2010), and the equations are written as follows:

$$\begin{array}{}\text{(12)}& {\displaystyle }{h}_{\mathrm{s},\text{eff}}^{i+\mathrm{1}/\mathrm{2},j}=max({\mathit{\eta }}_{\mathrm{s}}^{i,j},{\mathit{\eta }}_{\mathrm{s}}^{i+\mathrm{1},j})-max(z_{\mathrm{s}}^{i,j},z_{\mathrm{s}}^{i+\mathrm{1},j}),\text{(13)}& {\displaystyle }{h}_{\mathrm{s},\text{eff}}^{i,j+\mathrm{1}/\mathrm{2}}=max({\mathit{\eta }}_{\mathrm{s}}^{i,j},{\mathit{\eta }}_{\mathrm{s}}^{i,j+\mathrm{1}})-max(z_{\mathrm{s}}^{i,j},z_{\mathrm{s}}^{i,j+\mathrm{1}}),\end{array}$$

where η_s denotes water surface elevation (m), and ${\mathit{\eta }}_{\mathrm{s}}=h_{\mathrm{s}}+z_{\mathrm{s}}$. The friction terms S_fx and S_fy in Eq. (9) are evaluated with Manning's equation (Jain and Singh, 2005), and the equations of flow velocity at the boundary of Ω could be given as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}{u}^{i+\mathrm{1}/\mathrm{2},j}=& -\text{sgn}({\mathit{\eta }}_{\mathrm{s}}^{i+\mathrm{1},j}-{\mathit{\eta }}_{\mathrm{s}}^{i,j}){\displaystyle \frac{\mathrm{2}}{n_{\mathrm{s}}^{i+\mathrm{1},j}+n_{\mathrm{s}}^{i,j}}}\cdot \\ & (h_{\mathrm{s}}^{i+\mathrm{1}/\mathrm{2},j}{)}^{\mathrm{2}/\mathrm{3}}{\left|{\displaystyle \frac{{\mathit{\eta }}_{\mathrm{s}}^{i+\mathrm{1},j}-{\mathit{\eta }}_{\mathrm{s}}^{i,j}}{\mathrm{\Delta }x}}\right|}^{\mathrm{1}/\mathrm{2}},\end{array}\text{(15)}& {\displaystyle }\begin{array}{rl}{v}^{i,j+\mathrm{1}/\mathrm{2}}=& -\text{sgn}({\mathit{\eta }}_{\mathrm{s}}^{i,j+\mathrm{1}}-{\mathit{\eta }}_{\mathrm{s}}^{i,j}){\displaystyle \frac{\mathrm{2}}{n_{\mathrm{s}}^{i,j+\mathrm{1}}+n_{\mathrm{s}}^{i,j}}}\cdot \\ & (h_{\mathrm{s}}^{i,j+\mathrm{1}/\mathrm{2}}{)}^{\mathrm{2}/\mathrm{3}}{\left|{\displaystyle \frac{{\mathit{\eta }}_{\mathrm{s}}^{i,j+\mathrm{1}}-{\mathit{\eta }}_{\mathrm{s}}^{i,j}}{\mathrm{\Delta }x}}\right|}^{\mathrm{1}/\mathrm{2}},\end{array}\end{array}$$

where n_s is the surface roughness coefficient ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}$). Wall boundary conditions are applied, indicating that no surface runoff can flow through the simulation domain.

2.2.4 Subsurface runoff concentration

The subsurface runoff, which includes interflow and groundwater runoff, is routed with 2D linear reservoir equations. The governing equations of the original one-dimensional (1D) linear reservoir equations at the grid scale are given as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(\begin{array}{c}{O}_{\mathrm{i}}^{i,j}\\ O_{\mathrm{g}}^{i,j}\end{array}\right)=\left(\begin{array}{c}{Q}_{\mathrm{i},\text{up}}^{i,j}-Q_{\mathrm{i}}^{i,j}+{\mathit{\varphi }}_{\mathrm{i}}^{i,j}\\ Q_{\mathrm{g},\text{up}}^{i,j}-Q_{\mathrm{g}}^{i,j}+{\mathit{\varphi }}_{\mathrm{g}}^{i,j}\end{array}\right),\text{(17)}& {\displaystyle }\left(\begin{array}{c}{Q}_{\mathrm{i}}^{i,j}\\ Q_{\mathrm{g}}^{i,j}\end{array}\right)=\left(\begin{array}{c}-O_{\mathrm{i}}^{i,j}\ln \left(C_{\mathrm{i}}^{i,j}\right)/\mathrm{\Delta }T\\ -O_{\mathrm{g}}^{i,j}\ln \left(C_{\mathrm{g}}^{i,j}\right)/\mathrm{\Delta }T\end{array}\right),\end{array}$$

where O_i and O_g are the amount of interflow and groundwater storage, respectively (mm); Q_i and Q_g are the outflow intensity of the interflow and groundwater storage separately (mm s⁻¹); Q_i,up and Q_g,up are the summation of Q_i and Q_g from multiple upstream grids (mm s⁻¹); C_i and C_g are interflow and groundwater storage recession coefficients; ϕ_i and ϕ_g are source terms of interflow and groundwater storage, respectively (mm s⁻¹), which include the corresponding generated runoff and discharge into the channel; ${\mathit{\varphi }}_{\mathrm{i}}=R_{\mathrm{i}}-{\mathit{\epsilon }}_{\mathrm{i}}{Q}_{\mathrm{i}}$; and ${\mathit{\varphi }}_{\mathrm{g}}=R_{\mathrm{g}}-{\mathit{\epsilon }}_{\mathrm{g}}{Q}_{\mathrm{g}}$. When a channel is present in the grid, all outflows of the interflow and groundwater storage are assumed to drain into the channel, resulting in both ε_i and ε_g being equal to 1. In the no-channel case, ε_i and ε_g are set to 0.

To extend the linear reservoir equations from 1D to 2D, we follow the approach outlined by Liu et al. (2004). This involves decomposing the outflow intensity derived from the 1D method into x- and y-directional components based on the actual flow direction (see Fig. 2b). The state variables are then updated within the FVM framework. For clarity and simplicity, we will illustrate this extension using interflow given the similarity between the formulas for interflow and groundwater runoff concentration. The x- and y-directional components of Q_i, which maintain the total amount, could be expressed as follows:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& \left(\begin{array}{c}{Q}_{\text{i,x}}^{i,j}\\ Q_{\text{i,y}}^{i,j}\end{array}\right)=\\ & \left(\begin{array}{c}(\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{i}}^{i,j})Q_{\mathrm{i}}^{i,j}\sin {\mathit{\gamma }}^{i,j}/(\sin {\mathit{\gamma }}^{i,j}+\cos {\mathit{\gamma }}^{i,j})\\ (\mathrm{1}-{\mathit{\epsilon }}_{\mathrm{i}}^{i,j})Q_{\mathrm{i}}^{i,j}\cos {\mathit{\gamma }}^{i,j}/(\sin {\mathit{\gamma }}^{i,j}+\cos {\mathit{\gamma }}^{i,j})\end{array}\right),\end{array}\end{array}$$

where Q_i,x and Q_i,y are the x- and y-directional components of Q_i (mm s⁻¹); γ is the slope aspect, which is calculated clockwise, with 0 in the positive direction of the y axis; and $\mathrm{0}\mathit{°}\le \mathit{\gamma }<\mathrm{360}\mathit{°}$. The governing equation of the 2D linear reservoir method is given as follows:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}{O}_{\mathrm{i}}^{i,j}}{\mathrm{d}t}}=& -F_{\mathrm{i}}^{i+\mathrm{1}/\mathrm{2},j}+F_{\mathrm{i}}^{i-\mathrm{1}/\mathrm{2},j}-F_{\mathrm{i}}^{i,j+\mathrm{1}/\mathrm{2}}\\ & +F_{\mathrm{i}}^{i,j-\mathrm{1}/\mathrm{2}}+{\mathit{\varphi }}_{\mathrm{i}}^{i,j},\end{array}\end{array}$$

where $F_{\mathrm{i}}^{i+\mathrm{1}/\mathrm{2},j}$ and $F_{\mathrm{i}}^{i,j+\mathrm{1}/\mathrm{2}}$ denote the interflow intensity (mm s⁻¹) in the direction of the external normal vector at the eastern and northern boundaries of the control volume. The equations of $F_{\mathrm{i}}^{i+\mathrm{1}/\mathrm{2},j}$ and $F_{\mathrm{i}}^{i,j+\mathrm{1}/\mathrm{2}}$ could be expressed as follows:

$$\begin{array}{}\text{(20)}& \left(\begin{array}{c}{F}_{\mathrm{i}}^{i+\mathrm{1}/\mathrm{2},j}\\ F_{\mathrm{i}}^{i,j+\mathrm{1}/\mathrm{2}}\end{array}\right)=\left(\begin{array}{c}max(Q_{\text{i, x}}^{i,j},\mathrm{0})-min(Q_{\text{i, x}}^{i+\mathrm{1},j},\mathrm{0})\\ max(Q_{\text{i, y}}^{i,j},\mathrm{0})-min(Q_{\text{i, y}}^{i,j+\mathrm{1}},\mathrm{0})\end{array}\right).\end{array}$$

The same wall boundary conditions are applied to Eq. (19), which means that no subsurface runoff can flow through the simulation domain.

2.2.5 Channel concentration

The water within the channel is routed to the watershed outlet by using one-dimensional diffusion wave equations (Kazezyılmaz-Alhan and Medina, 2007), and the governing equations are as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial \mathit{U}}{\partial t}}+{\displaystyle \frac{\partial \mathit{F}\left(\mathit{U}\right)}{\partial x}}=\mathit{S}\left(\mathit{U}\right),\text{(22)}& {\displaystyle }\mathit{U}=\left[\begin{array}{c}{h}_{\mathrm{c}}\\ \mathrm{0}\end{array}\right]\mathit{F}=\left[\begin{array}{c}{h}_{\mathrm{c}}w\\ g{h}_{\mathrm{c}}^{\mathrm{2}}/\mathrm{2}\end{array}\right]\mathit{S}=\left[\begin{array}{c}{\mathit{\varphi }}_{\mathrm{c}}\\ g{h}_{\mathrm{c}}(S_{\text{oc}}-S_{\text{fc}})\end{array}\right],\end{array}$$

where h_c is channel water depth (m); w is channel flow velocity (m s⁻¹); S_oc is the channel bottom slope term (–); $S_{\text{oc}}=-\partial z_{\mathrm{c}}/\partial x$; z_c is the channel bottom elevation (m); S_fc is the channel friction term (–); and ϕ_c is the source term of the channel (m s⁻¹), including interflow recharge, groundwater recharge, and interaction discharge between the slope surface and channel.

The cell-centered FVM is used to discretize the 1D diffusion wave equations. The control volumes are constructed based on channel segment, with h_c being positioned at the center of the control volumes and with w being assigned along the boundaries. The mesh layout and stencils used are shown in Fig. 2c. The spatial discretization form of Eqs. (21) and (22) could be expressed as follows:

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{c}}^i}{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{1}}{A_{\mathrm{c}}^i}}[Q_{\text{c,up}}^i-Q_{\mathrm{c}}^i+{\mathit{\varphi }}_{\mathrm{c}}^i],\end{array}$$

where i is the index of the channel control volume (–), $A_{\mathrm{c}}^i$ is open-water surface area of the channel control volume (m²), $A_{\mathrm{c}}^i=\left[B\right(h_{\text{c,eff}}^{i-\mathrm{1}/\mathrm{2}})+B(h_{\text{c,eff}}^{i+\mathrm{1}/\mathrm{2}}\left)\right](\mathrm{\Delta }l{)}^i/\mathrm{2}$, B(⋅) is the formula of water surface width (m), (Δl)ⁱ denotes the length of the channel control volume (m), and $Q_{\mathrm{c}}^i$ is the channel discharge in the direction of the external normal vector at the eastern boundary of the channel control volume (m³ s⁻¹). Q_c,up is the summation of Q_c from multiple upstream channel segments. The formula of $Q_{\mathrm{c}}^i$ is derived based on Manning's equation, which is written as follows:

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{Q}_{\mathrm{c}}^i=& -\text{sgn}({\mathit{\eta }}_{\mathrm{c}}^{i+\mathrm{1}}-{\mathit{\eta }}_{\mathrm{c}}^i){\displaystyle \frac{\mathrm{2}}{n_{\mathrm{c}}^{i+\mathrm{1}}+n_{\mathrm{c}}^i}}\cdot \\ & {\displaystyle \frac{A(h_{\text{c,eff}}^{i+\mathrm{1}/\mathrm{2}}{)}^{\mathrm{5}/\mathrm{3}}}{\mathit{\chi }(h_{\text{c,eff}}^{i+\mathrm{1}/\mathrm{2}}{)}^{\mathrm{2}/\mathrm{3}}}}{\left|{\displaystyle \frac{{\mathit{\eta }}_{\mathrm{c}}^{i+\mathrm{1}}-{\mathit{\eta }}_{\mathrm{c}}^i}{(\mathit{\delta }l{)}^{i+\mathrm{1}/\mathrm{2}}}}\right|}^{\mathrm{1}/\mathrm{2}},\end{array}\end{array}$$

where A(⋅) denotes the formula of the cross-sectional area (m²), χ(⋅) is the formula of the wetted perimeter (m), η_c is the elevation of the channel water surface elevation (m), ${\mathit{\eta }}_{\mathrm{c}}=h_{\mathrm{c}}+z_{\mathrm{c}}$, n_c is the channel roughness coefficient ($\mathrm{s}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}$), and $(\mathit{\delta }l{)}^{i+\mathrm{1}/\mathrm{2}}$ is the distance between the center of the ith and i+1th channel control volume (m). The cross-section is generalized into a trapezoid, and the formulas of the cross-sectional hydraulic elements including A, B, and χ can be derived accordingly (see Eqs. S1, S2, and S3 in the Supplement). The effective water depth at the eastern boundary of the channel control volume h_c,eff is evaluated with the upwind scheme, which could be expressed as follows:

$$\begin{array}{}\text{(25)}& h_{\text{c,eff}}^{i+\mathrm{1}/\mathrm{2}}=\left\{\begin{array}{ll}{h}_{\mathrm{c}}^{i+\mathrm{1}}& {\mathit{\eta }}_{\mathrm{c}}^{i+\mathrm{1}}>{\mathit{\eta }}_{\mathrm{c}}^i\\ h_{\mathrm{c}}^i& {\mathit{\eta }}_{\mathrm{c}}^{i+\mathrm{1}}\le {\mathit{\eta }}_{\mathrm{c}}^i\end{array}\right..\end{array}$$

The exchange discharge between the slope surface and the channel is calculated based on their water surface elevation and Manning's equation (Shen and Phanikumar, 2010), which could be expressed as follows:

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{Q}_{\text{sc}}=& \text{sgn}(h_{\mathrm{s}}-h_{\text{c,eff}}){\displaystyle \frac{max(h_{\mathrm{s}},h_{\text{c,eff}}{)}^{\mathrm{5}/\mathrm{3}}\mathrm{\Delta }l}{n_{\mathrm{s}}}}\cdot \\ & {\left|{\displaystyle \frac{h_{\mathrm{s}}-h_{\text{c,eff}}}{\mathrm{\Delta }x/\mathrm{2}}}\right|}^{\mathrm{1}/\mathrm{2}},\end{array}\end{array}$$

where Q_sc is the exchange discharge between the slope surface and the channel (m³ s⁻¹); h_c,eff is the effective channel water depth (m), which exceeds the channel bank elevation z_bank (m); and $h_{\text{c,eff}}=max({\mathit{\eta }}_{\mathrm{c}}-z_{\text{bank}},\mathrm{0})$. The source term in Eq. (23) turns out to be

$$\begin{array}{}\text{(27)}& {\mathit{\varphi }}_{\mathrm{c}}=Q_{\text{sc}}+(Q_{\mathrm{i}}+Q_{\mathrm{g}})\mathrm{\Delta }{x}^{\mathrm{2}}/\mathrm{1000}.\end{array}$$

The upstream boundary of the headwater channel is the wall boundary condition, while the downstream boundary of the channel at the watershed outlet is the zero-depth gradient (ZDG) condition (Panday and Huyakorn, 2004). The equation of the ZDG condition is written as follows:

$$\begin{array}{}\text{(28)}& Q_{\text{ZDG}}={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{c}}^{\text{down}}}}{\displaystyle \frac{A(h_{\mathrm{c}}^{\text{down}}{)}^{\mathrm{5}/\mathrm{3}}}{\mathit{\chi }(h_c^{\text{down}}{)}^{\mathrm{2}/\mathrm{3}}}}{\left({\displaystyle \frac{\mathrm{2}{h}_c^{\text{down}}}{\mathrm{\Delta }{l}^{\text{down}}}}+S_{\text{oc}}^{\text{down}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where Q_ZDG is the outflow discharge at the watershed outlet (m³ s⁻¹), and the superscript “down” refers to the information stored at the channel control volume where the watershed outlet is located.

2.3 Numerical implementation

In the GXAJ model, the difference form equations for runoff generation from the original lumped XAJ model are directly applied at the grid scale. These equations are derived based on the time interval of input force, denoted as ΔT. However, owing to the Courant–Friedrichs–Lewy (CFL) stability criterion, the model time step Δt used for the runoff concentration module is usually smaller than ΔT. To address this discrepancy, a double-layer time loop strategy (DTLS) is adopted.

The DTLS consists of an outer-layer time loop that iterates forcing data and an inner-layer time loop that determines Δt for advancing the solution in time. Initially, the total simulation time T is read to define the duration of the outer-layer time loop. During this loop, forcing data are processed in chronological order. At the start of each outer-layer time loop, the forcing data and ΔT are loaded. The total amounts of three runoff sources (surface runoff, interflow, and groundwater runoff) are calculated using the difference form of the runoff generation equations. These amounts are then averaged to determine the intensities for three runoff sources based on ΔT. Following this, the inner-layer time loop begins, with ΔT serving as its time duration. The value of Δt is constrained by the CFL conditions and a user-defined maximum time step length (e.g., 10 min). The runoff concentration module is advanced in each inner time loop. To ensure that the end time of the last inner loop aligns with the start time of the next outer-layer iteration, Δt may be shortened as necessary. The solving process concludes once the total simulation time T is reached.

Given that the runoff generation and concentration processes are calculated separately and using equations in different forms, the numerical implementation is referred to as a loosely-coupled numerical implementation framework. The diagram of this framework is shown in Fig. 3a.

Figure 3 Diagram of the loosely coupled (a) and fully coupled (b) numerical implementation framework. t₁ and t₂ are two variables to store time. W_u, W_l, W_d, and S₀ are state variables related to the runoff generation process, while O_i, O_g, h_s, and h_c are state variables related to the runoff concentration process. The differences between the two frameworks are highlighted with dash boxes.

The mathematical equations of the TDD-XAJ model, after spatial discretization, form a set of ODEs. All state variables are assembled into a vector and are advanced simultaneously within a single time step to determine their values at the next time point. The equation of the coupled ODEs is written as follows:

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\mathrm{d}\mathit{Y}}{\mathrm{d}t}}=\mathit{F}(\mathit{Y},t),\end{array}$$

where Y is the vector of the model state variables; $\mathit{Y}={\left[\begin{array}{cccccccc}{W}_{\mathrm{u}}& W_{\mathrm{l}}& W_{\mathrm{d}}& S_{\mathrm{0}}& O_{\mathrm{i}}& O_{\mathrm{g}}& h_{\mathrm{s}}& h_{\mathrm{c}}\end{array}\right]}^{\top }$; and F is the function vector which represents the right-hand side of the ODEs, made up of a combination of Eqs. (1), (5), (11), (19), and (23). The first seven terms of Y are used to store the state variables of the slope, and the number of elements for each term corresponds to the total number of valid DEM grids. In contrast, the last term of Y stores the state variable related to the channel, and, thus, the number of elements for it matches the number of DEM grids occupied by river system.

This global coupling approach offers high numerical accuracy and flexibility for future expansion while ensuring the conservation and temporal continuity of all state variables (Qu and Duffy, 2007; Shu et al., 2020). Most importantly, it enables a fully coupled numerical implementation framework (see Fig. 3b). The DTLS is also used in this framework to address the mismatch between ΔT and Δt. The main differences between the fully and loosely coupled numerical implementation framework are as follows: (a) the total number of input forces – rather than the total number of runoff sources – is averaged in the outer time loop, and (b) both the runoff generation and concentration processes are calculated together using the differential form of equations.

The ODEs in both fully and loosely coupled numerical implementation frameworks are solved with the explicit Heun scheme (Clark and Kavetski, 2010). The combined restriction of CFL conditions for surface and channel flow is written as follows:

$$\begin{array}{}\text{(30)}& \mathrm{\Delta }{t}_{\text{max}}=\mathit{\alpha }min[\mathrm{\Delta }x/(\sqrt{(u^{\mathrm{2}}+v^{\mathrm{2}})}{)}_{\text{max}},\mathrm{\Delta }l/w_{\text{max}}],\end{array}$$

where Δt_max denotes the maximum time step without further user-defined or programmatic constraints (s), α is the CFL condition coefficient (–), and the rest subscript “max” denotes the maximum flow velocity across all corresponding elements.

2.4 Model parameters

The TDD-XAJ model has 15 tunable parameters, which is one more than in the original lumped Xinanjiang model (Zhao et al., 2023). The signification and value range of these parameters are listed in Table 1. Parameters that can be spatially uniform are consistent with those from the original lumped model. The methods for determining spatially distributed parameters are provided in Sect. S2 in the Supplement.

Table 1 Parameters and parameter ranges of the TDD-XAJ model.

3 Numerical experiments and model application

Two numerical experiments were conducted to evaluate the TDD-XAJ model, which was further applied in a typical humid watershed. The experiments compared 1D and 2D slope concentration methods, as well as loosely and fully coupled numerical implementation frameworks. These experiments were designed to demonstrate the theoretical effectiveness of the TDD-XAJ model, while the application was designed to assess the model's performance in real-world watersheds.

3.1 Numerical experiments

3.1.1 Slope concentration method comparison experiment

The performance of the 1D and 2D diffusion wave equations for surface runoff concentration and the performance of the 1D and 2D linear reservoir equations for subsurface runoff concentration were compared separately using two test cases. A total of eight simulations were conducted, categorized into surface and subsurface slope concentration comparison scenarios. The models used differ only in terms of their slope concentration methods. The evapotranspiration and runoff calculation module and the runoff separation module were disabled. The subsurface slope concentration module was turned off in the surface slope concentration comparison scenario and vice versa. The same 1D diffusion wave method was applied for channel concentration across all simulations.

The synthetic V catchment, first proposed by Overton and Brakensiek (1970), is commonly used to verify runoff concentration components in the hydrological modeling community (Kollet et al., 2017; Maxwell et al., 2014; Shen and Phanikumar, 2010; Shu et al., 2020). It consists of two symmetric hillslopes with a channel in between. Each hillslope has a length of 1000 m and a width of 800 m, with a surface roughness coefficient n_s of 0.015 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}$. The length and bottom slope of the channel S_oc are 1000 m and 0.02, respectively, and the channel roughness coefficient n_c is set to be 0.015 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}$. The channel has a uniform square cross-section, with a width measuring 20 m. This study used two test cases based on the synthetic V catchment, namely the single- and double-slope cases (Fig. 4), differentiated by the hillslope gradient in the y direction S_oy (parallel to the channel). Both cases have an x-directional gradient S_ox (perpendicular to the channel) of 0.05, while the S_oy is 0 for the single-slope case and 0.02 for the double-slope case.

Figure 4 Diagram of single-slope (a) and double-slope (b) synthetic-V-catchment test case.

In both test cases, precipitation lasted 180 min, with constant intensity during the first 90 min, followed by 90 min without precipitation, totaling 16.2 mm. Precipitation was uniformly distributed over the hillslopes, and the channel received no direct precipitation. Evapotranspiration was set to zero. Only the left-side hillslope was modeled considering the symmetry, and the hillslope was discretized into a structured grid with 5 m×5 m cells, totaling 32 000 cells, and the channel was discretized into 5 m long segments. The initial water depths of the hillslope and channel were set to zero. The wall boundary condition is applied to the hillslope, except for the boundary connected to the channel, allowing water to drain into the channel. The upstream boundary of the channel was set as the wall boundary condition, while the ZDG boundary condition was applied to the downstream outlet.

In the surface slope concentration comparison scenario, all precipitation was assumed to contribute to surface runoff with no losses. The simulated results using 1D and 2D diffusion wave equations were compared against analytical solutions and previous published results, respectively. The single-slope case has an analytical solution for hillslope and channel outflow, with a hillslope outflow duration of 180 min and a channel outflow duration of 90 min (Di Giammarco et al., 1996). No analytical solution exists for the double-slope case, and so the results of the Integrated Finite Difference Model (IFDM) were used as benchmarks (Di Giammarco et al., 1996). The Nash–Sutcliffe efficiency coefficient (NSE) (Nash and Sutcliffe, 1970) was used to assess the consistency of the simulation results with these benchmarks.

In the subsurface slope concentration comparison scenario, all precipitation was assumed to be transformed into interflow. The 1D and 2D linear reservoir equations are used to represent the interflow and groundwater runoff slope concentration processes, with interflow serving as the example in this experiment. The interflow storage recession coefficient C_i was set to 0.5, meaning that 50 % of the interflow reservoir storage would remain after an hour if no additional interflow entered.

3.1.2 Numerical implementation framework comparison experiment

The second numerical experiment aimed to investigate the effect of the numerical implementation framework on simulation results. It was conducted on the same synthetic-V-catchment test cases as the first numerical experiment. The runoff calculation, runoff separation, surface slope concentration, and channel concentration process are considered. Two different models emerged from the use of the loosely and fully coupled numerical implementation frameworks (see Sect. 2.3 for details), referred to hereinafter as the loosely and the fully coupled model.

The model parameters used in this experiment are presented in Table 1. A total of 500 sets of parameters were generated using the symmetric Latin hypercube (SLH) method (Gong et al., 2015), accounting for potential parameter variability. The parameters for hillslope and channel concentration were consistent with the previous experiments. The model parameters K_e, A_imp, K_i, K_g, C_i, and C_g were set to 0, ensuring that only surface runoff was involved, along with the derivation of the corresponding analytical solution. The total amounts of simulated surface runoff during the entire simulation using both loosely and fully coupled models were evaluated with mean absolute error (MAE) (Hodson, 2022) against the analytical value. The equation of this analytical value (Zhao et al., 2023) could be given as follows:

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{R}_{\mathrm{s}}^{\ast }=& P-{\displaystyle \frac{S_{\text{mm}}}{\mathrm{1}+\text{ex}}}[\mathrm{1}-(\mathrm{1}-P/S_{\text{mm}}{)}^{\mathrm{1}+\text{ex}}]\\ & -{\displaystyle \frac{W_{\text{mm}}}{\mathrm{1}+b}}[\mathrm{1}-(\mathrm{1}-P/W_{\text{mm}}{)}^{\mathrm{1}+b}]\\ & +{\displaystyle \frac{W_{\text{mm}}}{\mathrm{1}+\text{ex}+b}}[\mathrm{1}-(\mathrm{1}-P/W_{\text{mm}}{)}^{\mathrm{1}+\text{ex}+b}],\end{array}\end{array}$$

where $R_{\mathrm{s}}^{\ast }$ is the total surface runoff amount during the entire simulation (mm); P is the total amount of precipitation (mm); W_mm and S_mm are maximum single-point tension and free-water storage capacity, respectively (mm); $W_{\text{mm}}=W_{\mathrm{m}}(\mathrm{1}+b)/(\mathrm{1}-A_{\text{imp}})$; and $S_{\text{mm}}=S_{\mathrm{m}}(\mathrm{1}+\text{ex})$. The parameter S_m in each parameter set was adjusted to ensure that W_mm=S_mm, a precondition to be satisfied for the derivation of the analytical value of $R_{\mathrm{s}}^{\ast }$. The $R_{\mathrm{s}}^{\ast }$ for each parameter set is fixed at 16.2 mm for ease of comparison, and so the value of P associated with each parameter set is determined through a reverse calculation of Eq. (31).

However, there is no analytical solution for hillslope or channel outflow, making direct comparison challenging. To address this, we evaluated the convergence of the loosely coupled model by progressively reducing the time interval of input forcing ΔT. Theoretically, as ΔT decreases, the results of the loosely coupled model should converge to those of the fully coupled model. The initial ΔT was set to 90 min and reduced to 45 and 15 min for the loosely coupled model, while ΔT used for the fully coupled model was 90 min. The consistency of the loosely coupled model with the fully coupled model was evaluated using the MAE metric.

3.2 Model application

The TDD-XAJ model was applied in the Tunxi watershed (Fig. 5) to assess its performance. The Tunxi watershed is located at the headwater of the Xinanjiang River system in China, after which the XAJ model is named. The Tunxi watershed is a typical humid watershed, with annual precipitation of 1750 mm and an area of 2670 km². Its elevation ranges from 121 to 1614 m. The watershed is characterized by a hilly landscape comprising mountains, high and low hills, and intermountain basins. Due to the intense convective activity influenced by the unique terrain, precipitation can be heavy, resulting in rapid rises and falls due to floods.

Figure 5 Location and gauging station distribution of the Tunxi watershed (a) and the spatial discretization of the watershed, including channel and non-channel cells (b).

Daily-scale hydrological data were used in this study. The discharge and precipitation series were observed by the gauging network within the Tunxi watershed, which consists of the Tunxi hydrological station, the Yuetan hydrological station, and 23 precipitation stations (Fig. 5a). The duration ranges from 1 January 2007 to 31 December 2019, and the data source is the Annual Hydrological Report of China (Volume VI, Book XVIII). The daily pan evaporation series was collected from the daily surface climatological data for China (V3.0) within the same period. The spatial distribution of precipitation within the Tunxi watershed was obtained using the ordinary Kriging method provided by PyKrige (Murphy et al., 2024), while the watershed averaged value of pan evaporation was used for simulation.

A structured grid of 1000 m×1000 m was used for simulation (Fig. 5b). The digital elevation model (DEM) utilized was a 90 m resolution SRTMDEMUTM data product, while the land use information was obtained from the 30 m resolution GlobeLand30 (V2020) data product. To simplify the research and parameter calibration process, the parameters related to the runoff generation process were assumed to be spatially uniform and were calibrated manually. The surface roughness coefficient n_s was derived from the land use product (Shu et al., 2024). Based on DEM data and GIS analysis, 17 river channels were extracted for simulation. The channel roughness coefficient n_c was estimated based on the actual characteristics of the river channel. The range for n_c is between 0.025 and 0.040. The hydrometeorological data in the first year (1 January 2007 to 31 December 2007) were used for model spin-up, while the data in the following 7 years (1 January 2008 to 31 December 2014) were used for calibration, and the remaining data were used for validation. Model performance was evaluated with the NSE, Kling–Gupta efficiency coefficient (KGE), flood volume relative error (FVRE), and coefficient of determination (R²) (Jackson et al., 2019). The equation of the FVRE could be given as follows:

$$\begin{array}{}\text{(32)}& \begin{array}{rl}\mathrm{\text{FVRE}}=& \left({\sum }_{i=\mathrm{1}}^n{Q}_{\text{sim},i}-{\sum }_{i=\mathrm{1}}^n{Q}_{\text{obs},i}\right)/{\sum }_{i=\mathrm{1}}^n{Q}_{\text{obs},i}\\ & \cdot \mathrm{100}\mathit{\%},\end{array}\end{array}$$

where Q_sim,i and Q_obs,i represent simulated and observed discharge at time step i (m³ s⁻¹), and n is the length of the sequence.

4 Results and discussion

4.1 Slope concentration method comparison

The surface and subsurface slope concentration methods were compared for single- and double-slope cases, which resulted in four sets of simulation results. The simulated hydrographs of hillslope and channel outflow are presented in Fig. 6. The state variable distributions on the left hillside at the 60 min mark are shown in Fig. 7.

Figure 6 The hillslope and channel outflow hydrograph of the single- and double-slope synthetic-V-catchment test cases in the surface and subsurface comparison scenarios. In the surface comparison scenario (a–d), the one-dimensional (1D) and two-dimensional (2D) diffusion wave equations (DW) were evaluated, and the analytical solution for the single-slope case and the IFDM solution for the double-slope case were derived from Di Giammarco et al. (1996). The 1D and 2D linear reservoir equations (LR) were compared for the subsurface comparison scenario (e–h).

Figure 7 The spatial distribution of surface water depth (h_s) and interflow storage (O_i) on the left-side hillslope in both the single-slope (a, b) and double-slope (c, d) synthetic-V-catchment test cases at the 60 min mark. The state variable distributions shown are simulated using two-dimensional (2D) slope concentration methods. The corresponding results of 1D methods are identical to those obtained from the single-slope case simulated with 2D methods, regardless of the test case used. For a clear comparison, the spatial distribution of O_i in the upper-left corner has been zoomed.

As shown in Fig. 6a and b, the simulation results of the 1D and 2D diffusion wave methods coincide with each other, and this can be attributed to the consistency of the actual overland flow direction with that analyzed by the D8 method, which is perpendicular to the channel. The NSE values for the hillslope and channel outflow exceed 0.999 when compared with the analytical solution, demonstrating the applicability of 1D and 2D diffusion wave methods in the single-slope case. Significant differences exist in the hydrographs corresponding to the 1D and 2D diffusion wave methods in the double-slope case, with a notably pronounced lag in the channel outflow for the 1D diffusion wave method, as presented in Fig. 6d. It is worth noting that, for the 1D diffusion wave method, the flow directions analyzed using the D8 method in both single- and double-slope cases are identical, consistently with the simulation results. Using the IFDM solution as a comparative reference, the NSE values of the hillslope outflow are 0.991 and 0.998 for the 1D and 2D diffusion wave methods, while the NSE values for channel outflow are 0.865 and 0.992, respectively. The poor performance of the 1D diffusion wave method can be explained by its neglect of the y-directional flow component, which is rooted in its inability to capture the surface microtopography. For further illustration, the surface water depth of the left-side hillslope in the double-slope case for the 1D diffusion wave method is shown in Fig. 7a, and the corresponding result of the 2D diffusion wave method is displayed in Fig. 7c.

The subsurface slope concentration comparison scenario shows the same pattern; i.e., the simulation results of the 1D and 2D linear reservoir method are the same in the single-slope case, but there is a certain discrepancy in the double-slope case. In the single-slope case, the intensity of hillslope interflow outflow increases with the continuation of the precipitation in the first 90 min and remains constant in the following 90 min, while the channel outflow increased throughout the simulation, which is shown in Fig. 6e and f. The discrepancy in the double-slope case is also attributable to the mismatch between the flow direction used and the actual flow direction in the 1D linear reservoir approach. The outflow intensity from the hillslope and channel using the 2D linear reservoir method is slower than that of the 1D method, and the hillslope interflow outflow keeps decaying after the precipitation ceases, as depicted in Fig. 6g. The flow length of the subsurface flow path in the double-slope case is longer than that of the single-slope case, which means that the interflow slope concentration process is subject to a stronger storage effect, resulting in a slower outflow rate. Meanwhile, the simulation results of the 1D linear reservoir are identical in the single- and double-slope cases. Hence, the simulation results of the 2D linear reservoir method are more reasonable in the double-slope case. The comparison of interflow linear reservoir storage in the upper-left corner of the left-side hillslope at the 60 min mark also supports the ignorance of the y-directional interflow component for the 1D linear reservoir method in the double-slope case (Fig. 7b and d).

The result of this numerical experiment indicates that 1D slope concentration methods perform relatively poorly in the double-slope case, regardless of the surface or subsurface comparison scenario evaluated. This limitation arises partly from the inability of the D8 method to capture the actual flow direction accurately, highlighting the need for a more suitable single flow direction algorithm. Essentially, this reflects the shortcomings of the 1D method in characterizing complex terrains. In contrast, the 2D slope concentration methods exhibit good simulation accuracy in both the single- and double-slope cases. These 2D methods effectively represent flow direction by synthesizing flow velocities in both the x and y direction without reliance on specific flow direction algorithms. In summary, the findings from this numerical experiment underscore the advantages of using 2D slope concentration methods within the TDD-XAJ model.

4.2 Numerical implementation framework comparison

The loosely and fully coupled models were conducted on two V-catchment test cases with different ΔT values. Model parameters and the amount of precipitation were configured to ensure a total of 16.2 mm $R_{\mathrm{s}}^{\ast }$ was generated. The MAE metric was employed to evaluate the simulated model fluxes of the loosely coupled model against the analytical values or the solutions from the fully coupled model. A total of 500 parameter sets were utilized, resulting in 500 MAE values for each model flux, as detailed in Table 2.

Table 2 MAE statistics of model fluxes in numerical implementation comparison experiment.

^a The results of the fully coupled model are used as references to calculate the MAE values for hillslope and channel outflow, and so the corresponding value is empty.

As shown in Table 2, the maximum and average MAE values of $R_{\mathrm{s}}^{\ast }$ simulated by the fully coupled model are $\mathrm{4.19}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{2.84}\times {\mathrm{10}}^{-\mathrm{4}}$ mm, respectively, indicating a good alignment with the analytical solution. It is worth noting that the analytical solution of $R_{\mathrm{s}}^{\ast }$ using the same parameter set should be identical for both single- and double-slope cases. When ΔT is set to 90 min, the maximum and average MAE values of $R_{\mathrm{s}}^{\ast }$ simulated by the loosely coupled model are 5.93 and 4.57 mm, which account for 36.6 % and 28.2 % of their analytical value. As ΔT is reduced from 90 to 15 min for the loosely coupled model, the maximum and average MAE values for $R_{\mathrm{s}}^{\ast }$ become 2.0 % and 0.9 % of their analytical value, demonstrating a significant decrease trend but remaining relatively large compared to the fully coupled model.

As ΔT decreases, the MAE values for both hillslope and channel outflow as simulated by the loosely coupled model show a consistent downward tendency for both single- and double-slope cases. These MAE values, benchmarked against the fully coupled model's simulation results, suggest that the outputs of the loosely coupled model are gradually converging towards those of the fully coupled model. As an example, we illustrated the outflow hydrograph of various model configurations in the double-slope case (Fig. 8), using the first set of the SLH-generated parameters.

Figure 8 The simulated hillslope outflow (a) and channel outflow (b) of the loosely coupled (LC) model and the fully coupled (FC) model in the double-slope case. The time in parentheses denotes the time interval of the input forces (ΔT) used for the model. The FC model uses 90 min only, while the LC model uses 90, 45, and 15 min for the convergence test. The parameters used are the first set of the parameters sampled with the SLH method based on the ranges listed in Table 2. The inflection point at the 75 min mark in the hillslope outflow hydrograph of the LC model (using ΔT=15 min) is highlighted with a red-dashed box, along with time series of the first-order numerical derivative of hillslope outflow (ΔQ).

As shown in Fig. 8, the hillslope and channel outflow hydrograph of the loosely coupled model increasingly converge with those of the fully coupled model as ΔT decreases. Consequently, the impact of the numerical implementation framework on simulation results is analyzed based on the solution of the loosely coupled and fully coupled model (ΔT=90 min for both models). The average MAE values for channel outflow are 2.01 and 2.06 m³ s⁻¹ in single- and double-slope cases, respectively, which are lower than those for hillslope outflow. This could be explained by the cancel-out phenomenon of numerical error. The numerical error alternates between positive and negative values, and delayed effects are introduced into the runoff concentration process. The channel outflow at a specific moment is influenced by the cumulative effects of the numerical error from the current and previous time steps, resulting in a reduction in the magnitude of the numerical error. Although the average MAE value for channel outflow shows a certain decrease, its overall magnitude remains comparable to that of hillslope outflow, indicating that the numerical errors still require control. Additionally, the mean MAE values for hillslope outflow and channel outflow in the double-slope test case are 1.87 and 1.72 m³ s⁻¹, higher than those for the single-slope test case. This could be attributed to the concentration of overland flow on the downstream side of the channel in the double-slope case, where the overall numerical error is more pronounced. This suggests that more complex terrain may be more prone to amplifying numerical errors, a phenomenon that, to the best of our knowledge, has not been reported in previous studies.

It can also be seen from Fig. 8 that the simulation results of the loosely coupled model show a steady-state phenomenon when ΔT equals 90 and 45 min, and there is a clear inflection point in the simulated hillslope outflow hydrograph when ΔT equals 15 min. The false steady outflow phenomenon and inflection point observed indicate a potential limitation in capturing the transient behaviors of the loosely coupled numerical implementation framework and highlight the advantages of the fully coupled model's continuous updating mechanism.

This experiment highlights that significant numerical error arises from the loosely coupled numerical implementation framework previously employed. Additionally, the numerical error may be exacerbated owing to a potential terrain amplification effect. In contrast, the fully coupled numerical framework can effectively control numerical errors and capture transient behaviors. In summary, the results of this numerical experiment emphasize the benefits of utilizing purely differential mathematical equations and a fully coupled numerical framework within the TDD-XAJ model.

4.3 Model application result in Tunxi watershed

The TDD-XAJ model was applied in the Tunxi watershed, a typical watershed in a humid area, to simulate daily hydrography. The calibrated parameters of the model are listed in Table 3. The simulated result was evaluated with four performance metrics annually against the observations, which are listed in Table 4.

Table 3 Model parameters of the TDD-XAJ in the Tunxi watershed.

Table 4 Annually evaluated simulation performance metrics of the TDD-XAJ model in the Tunxi watershed.

For Tunxi station at the outlet of the Tunxi watershed, Table 4 indicates that the values of the FVRE metric are all within ±20 %, with the absolute values of the FVRE ($\left|\text{FVRE}\right|$) averaging 8.3 % and 4.5 % for the calibration and validation period, respectively. In terms of hydrograph evaluation, the average values of the NSE and KGE are 0.88 and 0.83 for the calibration period and 0.87 and 0.76 for the validation period; these are slightly better for the calibration period than for the validation period. The minimum value of R² is 0.83 for all years, and the average value for all years is 0.90. In a direct comparison, Tong (2022) conducted a similar daily simulation in the same watershed using the GXAJ model, reporting average NSE and $\left|\text{FVRE}\right|$ values of 0.85 % and 11.0 % between 2008 and 2017, respectively. In contrast, the TDD-XAJ model achieved average values of 0.87 for the NSE and 7.7 % for $\left|\text{FVRE}\right|$ in the same period. For Yuetan station within the Tunxi watershed, Table 4 shows that FVRE metric values remain within ±20 %. The average $\left|\text{FVRE}\right|$ is 7.3 % and 4.8 % for the calibration and validation periods, respectively. Meanwhile, the average value of the NSE is 0.82 for the calibration period and 0.84 for the validation period, and the average KGE is 0.80 and 0.77 for the calibration period and validation periods, respectively. The average R² across all years is 0.86. Figure 9 provides an example of the simulated hydrograph at the Tunxi and Yuetan stations of the TDD-XAJ model in 2008.

Figure 9 The simulated hydrograph at Tunxi station (a) and Yuetan station (b) of the Tunxi watershed in 2008 using the TDD-XAJ model.

We further analyzed the spatial simulation capability of the TDD-XAJ model using precipitation events from 6–8 June 2008, which occurred just before the flood peak and exhibited distinct spatial patterns (Fig. 10a–c). On 6 June, almost no rainfall was recorded; on 7 June, precipitation was concentrated in the upper-left corner; and on 8 June, it was focused in the lower-right corner. The simulation results for tension and free-water storage reflect these patterns. Specifically, on 7 June (Fig. 10d and g), the overall tension water storage was unsaturated, with minimal free-water storage. Subsequently, saturation of tension water occurred in the upper-left part of the watershed, and free-water storage increased significantly in areas with heavier rainfall (Fig. 10e and h). On 8 June, the tension water in the lower-right corner became fully saturated due to the coverage of precipitation (Fig. 10f). For free water (Fig. 10i), the storage in the same area also approached saturation, whereas the saturated areas in the upper part quickly dissipated due to the faster recession of free water.

Figure 10 The spatial distribution of precipitation for three days in 2008 (a–c), along with the spatial distribution of tension (d–f) and free-water storage (g–i) simulated by TDD-XAJ model for the three days.

The application results indicate a good and consistent agreement between the observation and the simulations of the TDD-XAJ model in the Tunxi watershed, demonstrating a slight improvement in performance metrics compared to the GXAJ model, particularly for flood volume. Additionally, the spatial analysis of simulation results demonstrates that the TDD-XAJ model captures spatial variability in water storage effectively in response to varying precipitation patterns.

4.4 Computational efficiency of the model

The implementation of the TDD-XAJ model is detailed in Appendix A. The Python language was chosen for its robust ecosystem in hydrology and Earth system sciences (Stacke and Hagemann, 2021; Murphy et al., 2024), which facilitates rapid development and experimentation. While Python, as an interpreted language, is generally less computationally efficient than compiled languages, it benefits from libraries like NumPy, which rely on compiled languages for underlying computations. To enhance computational performance, we leveraged NumPy's vectorized operations and optimized functions, which are generally faster than native Python code.

To compare the computational efficiency of the 1D and 2D slope runoff concentration methods, we repeated the models used in the corresponding comparison experiment 10 times and recorded the runtime for each. The average runtime values, obtained using an AMD 5950X CPU (3.4 GHz) in single-core mode, are summarized in Table 5.

Table 5 Average runtime recorded in the slope runoff concentration method comparison numerical experiment.

As shown in Table 5, replacing 1D slope runoff concentration methods with 2D methods incurs minimal computational overhead while potentially enhancing efficiency. For the diffusion wave methods used for surface runoff concentration, the runtime for the 2D form increases by 2.1 % and 5.2 % compared to the 1D form in the single-slope and double-slope cases, respectively. Conversely, for the linear reservoir method applied to subsurface runoff concentration, switching from 1D to 2D reduces runtime by at least 45.6 % across both test cases. The numerical computation procedure for the slope runoff concentration methods mainly involves two steps: flux calculation and a state variable update. In the flux calculation step, the 1D and 2D diffusion wave methods require calculation once per cell (based on GIS-derived flow direction) or twice (in both the x and y directions), respectively. For the linear reservoir method, both the 1D and 2D forms require calculation only once per cell, and the outflow flux evaluated is further decomposed into x- and y-directional components based on the slope aspect for the 2D form. When updating state variables, both the diffusion wave and linear reservoir methods in their 1D form require logical judgments based on GIS-derived flow direction to identify the inflow cells. In their 2D form, this can be directly determined using cell indexing. The element-wise logical judgment is relatively time-consuming and is comparable to a one-time flux calculation per cell. Therefore, for the diffusion wave method, transitioning from 1D to 2D form does not result in a significant increase in runtime, but, for the linear reservoir method, the computational overhead is effectively reduced.

4.5 Limitations and future research

The main limitation of the TDD-XAJ model is that it addresses only the numerical errors on the timescale from the ODE's numerical solution while neglecting errors arising from the spatial discretization of the PDE via the FVM method. Future research should consider using methods like manufactured solutions (Bisht and Riley, 2019) as an alternative for evaluating numerical errors arising from both spatial and temporal discretization when exact solutions are difficult to obtain. Additionally, exploring more advanced spatial discretization techniques, such as second-order FVM or the discontinuous Galerkin method (Shaw et al., 2021), along with more sophisticated temporal integration methods, could help control numerical errors more effectively.

In terms of computational efficiency, although we make extensive use of the NumPy library, which is generally faster than native Python code, we recognize that NumPy may still be slower in certain operations (e.g., element-wise calculations) compared to compiled languages. This is a well-known limitation within the Python community, and solutions like just-in-time (JIT) compilation have been proposed (Lam et al., 2015), which convert frequently executed script code into machine code with further automatic optimizations. Although this paper primarily focuses on presenting the theoretical aspects of the TDD-XAJ model, we plan to optimize the code, including the implementation of parallelization, in future work.

It is essential to validate the model with more real-world watershed data and to evaluate its uncertainty, particularly for watersheds with varying scales, diverse underlying physical conditions, and hydrological data at different temporal resolutions. The differential-form mathematical equations established for the TDD-XAJ model provide a solid foundation for future research, including combining deep learning for better model parameterization and process understanding (Höge et al., 2022; Li et al., 2024).

5 Summary and conclusions

The Xinanjiang model has been widely applied in China, and its modeling philosophy is referenced in the development of several other hydrological models commonly used worldwide. However, the development of a distributed version of the XAJ model has been somewhat stagnant. In this study, we propose a new grid-based distributed Xinanjiang model (TDD-XAJ). The features of the TDD-XAJ model include the usage of 2D slope concentration methods, the establishment of purely differential-form mathematical equations, and a fully coupled numerical framework for model solution. The performance of the TDD-XAJ model was tested and compared with the existing distributed Xinanjiang model (GXAJ, proposed by Yao et al., 2009), leading to the following conclusions:

1. The TDD-XAJ model replaces 1D diffusion wave equations for surface slope concentration with 2D diffusion wave equations; the new 2D linear reservoir equations were derived for subsurface slope concentration. The combination of the 2D methods forms the slope concentration module of the model. Comparisons conducted on two V-catchment test cases reveal that the 2D methods enhance the model's ability to capture the microtopography of complex terrain. Additionally, using 2D methods eliminates the need for careful selection of a single flow direction algorithm, which is required in 1D methods.

2. The ODE form equations describing the runoff generation process are integrated with 2D slope and 1D channel concentration equations. This integration establishes a purely differential-form mathematical framework of the TDD-XAJ model. The framework facilitates a fully coupled numerical implementation, resulting in a system of ODEs that can be solved simultaneously. A comparison experiment of the numerical implementation indicates that numerical error introduced by the previous loosely coupled manner in the GXAJ model is evident and likely to be amplified by complex terrains, while the use of the fully coupled manner in the TDD-XAJ model can effectively mitigate this issue.

3. The TDD-XAJ model was applied in the Tunxi watershed, a typical humid study area in China, and its results were compared with those from the GXAJ model. The findings indicated that the TDD-XAJ model outperforms the GXAJ model in terms of key hydrological skill metrics, particularly in its improved ability to simulate flood volume. These results preliminarily demonstrate the effectiveness and applicability of the TDD-XAJ model for real-world watershed simulations.

Appendix A: Implementation details of the TDD-XAJ model

The TDD-XAJ model is implemented in Python 3.9, taking advantage of its extensive ecosystem, as well as its cross-platform capabilities. To maintain readability and facilitate understanding, only essential libraries are used, ensuring that the model is both accessible and easily extendable in the future. A list of these libraries is provided in Table A1. NumPy plays a crucial role, providing vectorized operations and optimized computational functions that enhance performance compared to native Python code. Other libraries serve supporting roles, such as model configuration and handling file input and output.

Table A1 Libraries utilized in the implementation of the TDD-XAJ model.

Appendix B: Inputs and outputs of the TDD-XAJ model

The inputs and outputs of the TDD-XAJ model are listed in Table B1. Its inputs consist of four categories: (1) model configuration dictionary, (2) spatiotemporal multidimensional arrays of precipitation and evaporation as forcing inputs, (3) raster datasets providing spatial attributes (including elevation, river network, aspect, etc.), and (4) lookup tables storing hydraulic parameters for surface runoff and channel concentration. The outputs comprise spatiotemporal distributions of state variables and model fluxes at pre-defined temporal resolutions (instantaneous or time-averaged values). Compared to the original XAJ model, which is a lumped hydrological model, its inputs are watershed-scale parameters and time series for precipitation and evaporation, with output mainly including the simulated discharge hydrograph. The TDD-XAJ model shares largely similar input and output information with the existing distributed XAJ model; the aspect data introduced by the 2D linear reservoir method can be derived from the elevation raster.

Table B1 Detailed information on input and output files of the TDD-XAJ model.