# A 2D shallow water flow model with 1D internal boundary condition for subgrid-scale topography

Pablo Vallés, J. Fernández-Pato, Mario Morales-Hernández, I. Echeverribar, Pilar García-Navarro

► To cite this version:

Pablo Vallés, J. Fernández-Pato, Mario Morales-Hernández, I. Echeverribar, Pilar García-Navarro. A 2D shallow water flow model with 1D internal boundary condition for subgrid-scale topography. Advances in Water Resources, 2024, 189, pp.104716. 10.1016/j.advwatres.2024.104716. hal-04947036

HAL Id: hal-04947036

https://hal.science/hal-04947036v1

Submitted on 6 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# A 2D shallow water flow model with 1D internal boundary condition for subgrid-scale topography

P. Vallés<sup>a,*</sup>, J. Fernández-Pato<sup>b</sup>, M. Morales-Hernández<sup>a</sup>, I. Echeverribar<sup>a,c</sup>, P. García-Navarro<sup>a</sup>

## ARTICLE INFO

### Keywords:

- Subgrid model

- Internal boundary condition

- Numerical simulation

- Finite volume

- Shallow water equations

## ABSTRACT

In this work, a dynamic internal boundary condition is used as subgrid model in a two-dimensional (2D) model based on the shallow water equations in order to model narrow regions in the domain. In this way, computational savings are sought, since it is not necessary to discretize these regions with cells of reduced size. The new internal boundary condition simplifies other works where 1D-2D coupled models were presented, since the 1D model is a subgrid for the 2D mesh, so the coupling between both models is simple and direct. The coupling is performed using mass conservation, simplifying the calculation in the transfer between both models. Test cases are studied to validate the implemented boundary condition, and a mountain catchment as a realistic case. The results obtained with a fully 2D mesh and a 2D mesh with rills in narrow regions are very similar, with a large reduction in computational cost when using rills, both in test cases and in the realistic case. Thus, the use of the implemented internal boundary condition is an effective tool to study regions with narrow regions by reducing the computational cost with little loss of accuracy in the results.

## 1. Introduction

The reduction of computational cost without loss of accuracy in the numerical results when simulating shallow water free surface flows is one of the most important challenges of computational hydraulics. In particular, when trying to reproduce events on a large spatial domain the efficiency of the method becomes crucial to make a computational tool affordable and with a practical use. The 2D shallow water system of equations has been considered suitable to reproduce the flow behavior in many situations (Sanders et al., 2010; Masoero et al., 2013; Defina, 2000; Costabile et al., 2017; Vacondio et al., 2016; Briganti and Dodd, 2009). The increasing availability of digital topographic data in recent years provides this type of models with a wider scope of application. However, this model may involve an excessive computational cost, specially due to the necessity of small cells to achieve the required resolution of bed topography (Bomers et al., 2019; Hu et al., 2019). This drawback of 2D models has caused during years an extended use of 1D models (Masood and Takeuchi, 2012; Aggett and Wilson, 2009; Petaccia et al., 2012) that appear advantageous when simulating complex channel nets. Despite the high computational efficiency, the use of 1D models often implies a reduced accuracy in the results compared to those provided by 2D models (Horritt and Bates, 2002; Costabile et al., 2015; Betsholtz and Nordlöf, 2017). For this reason, some researchers

have proposed 1D-2D coupled models to combine 2D representation on the main part of the domain and 1D numerical schematization for the main channels (Verwey, 2001; Dhondia and Stelling, 2002; Marin and Monnier, 2009; Gejadze and Monnier, 2007), leading not only to a reduction of the simulation time, but also to a more accurate representation of the river channel. Although this approach has been successfully used for many flood studies, it is generally time-consuming in setting up the initial model and the accuracy of predictions varies with the way in which the floodplain is discretized.

Some authors (Lin et al., 2006; Villanueva and Wright, 2006) suggest using only a 1D model for predicting flow velocity and water level in the main river network. However, if a breach in a river embankment causes extensive flooding, the flow characteristics would no longer remain strictly 1D. In such scenarios, a 2D model is employed to estimate flow velocity and inundation levels in the flooded area. These models are interconnected using a weir equation, which determines the flow volume from the 1D domain to the 2D domain based on the water level difference.

Another approach to coupling 1D-2D hydrodynamic models involves transforming 2D quantities into 1D equivalents by averaging the terms along cross-sections and enforcing continuity at the interfaces. Subsequently, an iterative procedure is undertaken within subdomains

<footnote><sup>a</sup> I3A, University of Zaragoza, Zaragoza, Spain</footnote>

<footnote><sup>b</sup> CSIC - Estación Experimental de Aula Dei (EEAD), Zaragoza, Spain</footnote>

<footnote><sup>c</sup> Hydronia Europe S.L., Madrid, Spain</footnote>

<footnote>* Corresponding author.
E-mail address: pvalles@unizar.es (P. Vallés).</footnote>
to solve the coupled 1D–2D problem (Miglio et al., 2005). This technique proves to be a reliable strategy when a suitable subdomain is selected, typically for simpler configurations such as straight channels or river bifurcations.

Recent studies propose more advanced methods to integrate both models. For instance, they establish internal connections by coupling the 1D node with the center of the 2D grid cell (Li et al., 2009), considering numerical fluxes from each model, and introducing corrections in the transfer of momentum quantity due to swirls (Finaud-Guyot et al., 2011).

Most of these coupled model approaches, developed from existing 1D and 2D models, require a comprehensive examination of how each model perceives the coupling independently. The boundary conditions in each model play a crucial role in the modeling process, particularly because the end of the 2D domain interacts continuously with the 1D model, needing a careful consideration of boundary treatment. Moreover, other coupling strategies based on a mass conservation and a complete mass and momentum conservation are proposed in recent works in order to cover all possible flow situations and to approximate faithfully the results given by a fully 2D model (Morales-Hernández et al., 2013).

However, in general, the coupling between the 1D and 2D models is an added complexity when creating the mesh (Marin and Monnier, 2009; Morales-Hernández et al., 2013; Brunner, 2014). And, depending on the case studied, such coupling is difficult to achieve, as in the case of steeply sloping basins, where the one-dimensional nature of the flow remains in narrow regions of the domain such as gullies. For that reason, in the present work a 2D model is proposed that, through the use of a simplified one-dimensional internal boundary condition (IBC), models these regions in a subgrid style, saving a large number of two-dimensional cells, which would be necessary if these regions were fully discretized in the 2D model. Compared to previous works, this is a simpler strategy, since the 1D model is a subgrid for the 2D mesh, and therefore one-dimensional cells do not exist, thus avoiding coupling problems between both models. Subgrid-scale modeling (SGS) refers to the representation of small-scale physical processes occurring at length scales that cannot be adequately resolved on a computational mesh or which resolution involves high computational consumption (Li et al., 2011; Kitts et al., 2020), as in our case. The use of subgrid models in certain regions is justified by obtaining accurate data regardless of the 2D mesh resolution, since finer 2D meshes imply a higher computational cost (Sehili et al., 2014; Casulli, 2009). Unlike the 1D model proposed as the subgrid model in this work, there are other subgrid models that model the behavior of certain variables, such as porosity, by introducing an additional source term in certain cells of the mesh (Bates, 2000; Nihei et al., 2024). In this work, the coupling of both models is based on the conservation of mass, so the transfer between the 2D model and the simplified 1D model is significantly simplified when compared to other strategies that also consider the conservation of linear momentum. Both models are discretized using a conservative upwind cell-centered finite volume scheme based on Roe Riemann solver across the edges (Roe, 1981). The topography is represented only with DTM (Digital Terrain Model) in a triangular unstructured grid for the 2D model that contains the rill (simplified channel) as an internal boundary condition. Internal boundary conditions are frequently used to model certain hydraulic structures in rivers, due to the significant reduction in computational consumption without losing accuracy in the results (Echeverriar et al., 2019; Dazzi et al., 2020). In the present work, test cases are presented to evaluate the relative performance of the simplified 1D–2D model and the fully 2D in steady and unsteady cases with the focus on relative performance and computational saving. Finally, it is applied to a rainfall runoff case mountain catchment with complex topography where the numerical results of the proposed model are compared with those obtained with a fully 2D modelization. The technique proposed in this work aims to address the problem of the mesh refinement by drastically reducing

the number of cells in areas where (1) the flow is predominantly unidirectional and (2) the density of cells is necessarily high to be able to faithfully discretize the topography of the basin, such as the main channel of the basin or the most significant ravines. In order to do this, these regions will be modeled using 1D rills connected to the rest of the 2D domain, which will be discretized using an unstructured triangular mesh, similar to the one used for reference simulations.

## 2. Governing equations

The mathematical models in our case are the shallow water equations.

### 2.1. 2D shallow water equations

In search for a compromise between computational efficiency and accuracy, 3D models are not usual in overland flow, since the information provided by a depth averaged 2D model is sufficiently accurate for the vast majority of problems studied (Vreugdenhil, 1994), assuming that the vertical acceleration is not relevant. The resulting differential continuity equation or mass conservation equation relates the water depth h to the average velocity components u and v in the $\hat{x}$ and $\hat{y}$ directions, respectively:

$$ \frac{\partial h}{\partial t} + \frac{\partial(hu)}{\partial x} + \frac{\partial(hv)}{\partial y} = q_{\text{out}} - q_{\text{in}} + R - f \quad (1) $$

being the discharge per unit width in the $\hat{x}$ component, $q_x = hu$, and the discharge per unit width in the $\hat{y}$ component, $q_y = hv$. In addition, $q_{\text{out}}$ and $q_{\text{in}}$ are the outflow and inflow discharges to the rill per unit area, acting as volume sink or source of the 2D model, R is the rainfall intensity and f the soil infiltration rate, modeled by means of the Horton infiltration model (Horton, 1933). Horton model predicts the soil infiltration capacity $f_p$ by assuming an exponential decay with time, as follows:

$$ f_p = f_c + (f_0 - f_c)e^{-kt} \quad (2) $$

being $f_0$ and $f_c$ the initial and equilibrium infiltration capacities, respectively, and k a constant representing the rate of decrease in $f_p$.

Eq. (2) is valid as long as the surface remains ponded ($h > 0$). In this case, the actual infiltration rate is equal to the infiltration capacity ($f = f_p$). On the other hand, when rainfall is falling on soil with no surface water, an additional consideration should be made. If the rainfall intensity is lower than the soil infiltration capacity at a given time ($R \le f_p$), all the incoming rainfall infiltrates and thence the actual infiltration rate becomes limited by the rainfall intensity ($f = R$). If rainfall exceeds the soil infiltration capacity ($R > f_p$), the soil becomes ponded and Eq. (2) predicts the actual infiltration rate ($f = f_p$).

The differential equations of momentum incorporate the influence of the driving force of gravity and the force of friction:

$$ \frac{\partial(hu)}{\partial t} + \frac{\partial(hu^2)}{\partial x} + gh\frac{\partial h}{\partial x} + \frac{\partial(huv)}{\partial y} = gh(S_{ox} - S_{fx}) \quad (3) $$

$$ \frac{\partial(hv)}{\partial t} + \frac{\partial(huv)}{\partial x} + gh\frac{\partial h}{\partial y} + \frac{\partial(hv^2)}{\partial y} = gh(S_{oy} - S_{fy}) \quad (4) $$

where $S_{ox}$ and $S_{oy}$ are the components of the gradient of the bottom surface $z_b(x, y)$:

$$ S_{ox} = -\frac{\partial z_b}{\partial x}, \quad S_{oy} = -\frac{\partial z_b}{\partial y} \quad (5) $$

while $S_{fx}$ and $S_{fy}$ are the frictional slopes of an energy surface:

$$ S_{fx} = \frac{n^2 u \sqrt{u^2 + v^2}}{h^{4/3}}, \quad S_{fy} = \frac{n^2 v \sqrt{u^2 + v^2}}{h^{4/3}} \quad (6) $$

with n the Manning's coefficient, which indicates the resistance that a flow suffers due to the bed roughness, whose empirical values are provided in tabulated references and obtained experimentally (Bcc Waterways Program Staff et al., 2000; Arcement and Schneider, 1984).
<figure><img src="image_4.png" /><figcaption>Fig. 1. Main 1D flow variables (a) and the representation of the wetted cross-section area and the wetted perimeter (b).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Representation of the relevant variables in the 2D numerical scheme. (a) Detail of the triangular mesh focusing on cell i and its neighbors. (b) Sketch of the numerical contributions defined at the edges of cell i and their sense of propagation.</figcaption></figure>

The expressions (1), (3) and (4) in conservative form can also be written in matrix form:

$$ \frac{\partial U}{\partial t} + \frac{\partial F(U)}{\partial x} + \frac{\partial G(U)}{\partial y} = S(U) \qquad (7) $$

with:

$$ U = \begin{pmatrix} h \\ hu \\ hv \end{pmatrix}, \quad F(U) = \begin{pmatrix} hu \\ hu^2 + \frac{1}{2}gh^2 \\ huv \end{pmatrix}, \quad G(U) = \begin{pmatrix} hv \\ huv \\ hv^2 + \frac{1}{2}gh^2 \end{pmatrix}, $$

$$ S(U) = \begin{pmatrix} q_{\text{out}} - q_{\text{in}} + R - f \\ gh[S_{\text{ox}} - S_{\text{fx}}] \\ gh[S_{\text{oy}} - S_{\text{fy}}] \end{pmatrix} \qquad (8) $$

where U is the vector of conserved variables; F(U) and G(U) are the fluxes of U in $\hat{x}$ and $\hat{y}$, respectively; and S(U) represents the vector of source terms.

## 2.2. 1D shallow water equations

Disregarding the transverse velocity component, the equations for 1D flow in prismatic channels can be written in the form:

$$ \frac{\partial U}{\partial t} + \frac{\partial F(U)}{\partial s} = S(U) \qquad (9) $$

where s is the longitudinal coordinate of the 1D channel and:

$$ U = \begin{pmatrix} h \\ q \end{pmatrix}, \quad F(U) = \begin{pmatrix} q \\ \frac{q^2}{h} + \frac{1}{2}gh^2 \end{pmatrix}, \quad S(U) = \begin{pmatrix} q_{\text{in}} - q_{\text{out}} \\ g[h(S_0 - S_f)] \end{pmatrix} \qquad (10) $$

where q is the discharge per unit area and h is the water depth. Fig. 1a shows a sketch of a rectangular cross section with the involved distances including $H = h + z_b$. The discharges $q_{\text{in}}$ and $q_{\text{out}}$ are the discharges per unit area as in (8). They represent the volume exchange between the 1D model and the 2D model. $S_0$ is the slope of the bottom in the longitudinal direction of the channel:

$$ S_0 = -\frac{\partial z_b}{\partial s} \qquad (11) $$

and $S_f$ is friction slope, based here again on Manning's law:

$$ S_f = \frac{q|q|n^2}{h^2 R^{4/3}} \qquad (12) $$

where R is the hydraulic radius, which is defined as the ratio of the wetted cross-section area $A = hB$ to the wetted perimeter $P = 2h + B$ (see Fig. 1b) and n is the Manning's coefficient (Bcc Waterways Program Staff et al., 2000; Arcement and Schneider, 1984). The variables A and R are obtained considering that the channel has a constant width B, facilitating the calculation of the simplified 1D hydrodynamic model.

## 3. Numerical resolution

The technique used to solve systems of Eqs. (7) and (9) is an explicit upwind finite volume scheme, based on the Riemann Solver of Roe (Toro, 1997). This numerical scheme starts from (7) or (9) which in general form can be expressed as:

$$ \frac{\partial U}{\partial t} + \vec{\nabla} \cdot E(U) = S(U) \qquad (13) $$
<figure><img src="image_6.png" /><figcaption>Fig. 3. Representation of the 1D numerical scheme for a subcritical flow.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 4. 2D mesh with a rill as IBC (a) and zoom of the grid cells around the rill (b).</figcaption></figure>

where $E(U) = (F(U), G(U))^{\top}$ in a 2D model, where $F(U)$ and $G(U)$ are given by (8); and $E(U) = F(U)$ in a 1D model, with $F(U)$ given by the expression (10). Integrating (13) into a volume or cell $\Omega$ and applying the divergence theorem:

$$ \frac{d}{dt} \int_{\Omega} \mathbf{U} \, d\Omega + \oint_{\partial\Omega} \mathbf{E}(U) \cdot \hat{\mathbf{n}} \, dl = \int_{\Omega} S(U) \, d\Omega \quad (14) $$

where $\partial\Omega$ is the contour of the volume and $\hat{\mathbf{n}}$ is the unit vector normal to the volume $\Omega$ defined in the outgoing direction of the contour. The discrete formulation of (14) is then outlined.

### 3.1. 2D numerical scheme

The discrete formulation of (14) in a 2D domain can be written as:

$$ \Omega_i \frac{\mathbf{U}_i^{n+1} - \mathbf{U}_i^n}{\Delta t_{2D}} + \sum_{k=1}^{NE} (\delta\mathbf{E})_k \cdot \hat{\mathbf{n}}_k l_k = \sum_{k=1}^{NE} \mathbf{S}_k \quad (15) $$

where $\Omega_i$ is the area of the cell $i$, $n$ is the index associated to the current time state and NE is the number of neighboring cells to cell $i$. Considering triangular cells, NE = 3. Furthermore, $\delta\mathbf{E} = \mathbf{E}_j - \mathbf{E}_i$, where $\mathbf{E}_j$ is the flux value at cell $\Omega_j$, which shares wall $k$, of length $l_k$, with cell $\Omega_i$ with flux value $\mathbf{E}_i$. In order to clarify the variables just presented, Figs. 2a and 2b are shown.

Taking into account the characteristics of the system (13), the Jacobian matrix of the flow can be defined:

$$ \mathbf{J}_{n,2D} = \frac{\partial(\mathbf{E}(U) \cdot \hat{\mathbf{n}})}{\partial U} = \frac{\partial F(U)}{\partial U} n_x + \frac{\partial G(U)}{\partial U} n_y \quad (16) $$

which can also be defined locally and made diagonal:

$$ \bar{\mathbf{J}}_{n,2D} = \bar{\mathbf{P}}_k \bar{\Lambda}_k \bar{\mathbf{P}}_k^{-1} \quad (17) $$

where the variables with tilde represent the average value at each cell and where $\bar{\Lambda}_k$ is the diagonal matrix whose non-zero elements are the eigenvalues of the system, $\bar{\lambda}^m$,

$$ \bar{\lambda}_1 = \bar{u} \cdot \hat{n} - \bar{c}, \quad \bar{\lambda}_2 = \bar{u} \cdot \hat{n}, \quad \bar{\lambda}_3 = \bar{u} \cdot \hat{n} + \bar{c} \quad (18) $$

being $\tilde{u} \cdot \hat{n} = \tilde{u} n_x + \tilde{v} n_y$ and $\tilde{c}$ the celerity of the infinitesimal surface deformation waves, and $\bar{\mathbf{P}}$ is the matrix containing the eigenvectors of the system, $\bar{\mathbf{e}}^m$:

$$ \bar{\mathbf{e}}_1 = \begin{pmatrix} 1 \\ \tilde{u} - \tilde{c} n_x \\ \tilde{v} - \tilde{c} n_y \end{pmatrix}, \quad \bar{\mathbf{e}}_2 = \begin{pmatrix} 0 \\ -\tilde{c} n_y \\ \tilde{c} n_x \end{pmatrix}, \quad \bar{\mathbf{e}}_3 = \begin{pmatrix} 1 \\ \tilde{u} + \tilde{c} n_x \\ \tilde{v} + \tilde{c} n_y \end{pmatrix} \quad (19) $$

Starting from (15), the updating expression of the conserved variables at each cell $i$, taking into account the contributions $m$ for each of its walls $k$ is obtained (see Fig. 2b):

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t_{2D}}{\Omega_i} \sum_{k=1}^{3} \sum_{m=1}^{3} [(\tilde{\lambda}^m - \tilde{\gamma} \bar{\mathbf{e}}_k)^m l_k]^n \quad (20) $$

The time step $\Delta t_{2D}$ is calculated dynamically according to:

$$ \Delta t_{2D} = \text{CFL min}_{k,m} \left( \frac{\delta x_k}{\bar{\lambda}_k^m} \right) \quad (21) $$

where CFL is the Courant-Friedrichs-Lewy number, with $0 < \text{CFL} \le 1$, which guarantees stability in the numerical scheme (LeVeque, 1992), and where:

$$ \delta x_k = \min(\chi_i, \chi_j) \quad (22) $$

being:

$$ \chi_i = \frac{\Omega_i}{\max_k l_k} \quad (23) $$

According to (21), the time step depends on the dynamics of the problem through $\bar{\lambda}_k^m$ and on the cell size. Therefore, increasing the refinement of the computational mesh will imply a smaller time step, and therefore a higher computational cost (Hu et al., 2019). It is essential to optimize the mesh using internal boundary conditions (Echeverriar et al., 2019) that represent correctly the topography of each region, thus avoiding an unnecessary increase in computational cost.
<figure><img src="image_8.png" /><figcaption>Fig. 5. Real representation (a) and schematic representation (b) of 2D and 1D models in a generic case.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Representation of 2D cells and 1D cell.</figcaption></figure>

## 3.2. 1D numerical scheme

A Jacobian matrix can also be written for the 1D model:

$$J_{n,1D} = \frac{\partial F(U)}{\partial U} = \begin{pmatrix} 0 & 1 \\ c^2 - u^2 & 2u \end{pmatrix} \quad (24)$$

whose eigenvalues are:

$$\lambda^{-} = u - c, \quad \lambda^{+} = u + c \quad (25)$$

with u and c given by:

$$u = \frac{Q}{Bh}, \quad c = \sqrt{\frac{A}{B}} = \sqrt{gh} \quad (26)$$

where A is the wetted cross-section area and B is the width of the free surface (see Fig. 1a).

The 1D numerical scheme is formulated for the update of a cell i between time t^n and time t^{n+1}, where the right and left contributions of fluxes and source terms ȳ of the cell interfaces are taken into account (see Fig. 3):

$$U_i^{n+1} = U_i^n - \frac{\Delta t_{1D}}{\Delta s_i} \left[ \sum_{m=1}^{2} (\tilde{\lambda}^+ \tilde{\gamma} \tilde{e})_{i-1/2}^m + \sum_{m=1}^{2} (\tilde{\lambda}^- \tilde{\gamma} \tilde{e})_{i+1/2}^m \right]^n \quad (27)$$

where ȳ and ȳ̃ are respectively the local averaged eigenvalues and eigenvectors of the Jacobian matrix of the flow evaluated on the wall, and where Δs_i is the i cell size (see Fig. 3). The eigenvalues are separated into their positive and negative parts:

$$\tilde{\lambda}_{i+1/2}^{\pm m} = \frac{1}{2}(\tilde{\lambda} \pm |\tilde{\lambda}|)_{i+1/2}^{m} \quad (28)$$

In (27), numerical stability requires that the time step Δt_1D is constrained by the CFL condition:

$$\Delta t_{1D} = \text{CFL min}_{m,k,l} \left( \frac{\Delta s_i}{|\tilde{\lambda}_k^m|} \right) \quad (29)$$

where 0 < CFL ≤ 1 (LeVeque, 1992).

## 4. Rill internal boundary condition

As in any system of evolutionary differential equations, it is necessary to consider initial conditions and boundary conditions. Within the boundary conditions, there are the so-called internal boundary conditions, which impose certain restrictions at cells in the interior of the domain (Echeverribar et al., 2019; Dazzi et al., 2020). The rill internal boundary condition will control the exchange of volume between the 2D domain and the simplified 1D model.

### 4.1. 2D mesh and 1D sub-grid

Unlike previous works (Marin and Monnier, 2009; Morales-Hernández et al., 2013; Brunner, 2014), the way to introduce the rill internal boundary condition is not by coupling 1D and 2D meshes, but to create a subgrid in the narrow region where the rill is located (see Fig. 4a). In this way, the 1D cells do not physically exist, but form a subgrid where the size of the 1D cells come according to the size of the pairs of 2D cells. There is no constraint but, for simplicity, in the examples shown in the present work, the size of the 1D cell has been chosen equal to that of the 2D cells side (see Fig. 4b). Therefore, the coupling of both models is greatly simplified, being always a left and right lateral coupling between the 2D cells of each pair and the 1D cell that is located between them.

### 4.2. Transferred volume in rill internal boundary condition

The volume transfer from the 2D regions to the one-dimensional region is based on the flow surface levels H = h + z_b at the cells on the right (H_{2D, right}) and left (H_{2D, left}) sides of the rill and within the rill
<figure><img src="image_10.png" /><figcaption>Fig. 7. Representation of the 1D numerical scheme for a supercritical flow.</figcaption></figure>

itself ($H_{1D}$). Thus, the following expression can be established for the unit discharge into the rill depending on whether the rill side, treated as a weir, is submerged or not (Echeverribar et al., 2019; Sotelo, 1977):

$$ q_{\text{in, left}} = \begin{cases} 0 & \text{if } H_{\text{w, left 2D}} \le 0, H_{\text{w, left 1D}} \le 0 \\ \frac{2}{3} \sqrt{2g} C_d H_{\text{w, left 2D}}^{3/2} & \text{if } H_{\text{w, left 2D}} > 0, H_{\text{w, left 1D}} \le 0 \\ \frac{2}{3} \sqrt{2g} C_d H_{\text{w, left 2D}}^{3/2} \left(1.0 - \left(\frac{H_{\text{w, left 1D}}}{H_{\text{w, left 2D}}}\right)^{3/2}\right)^{0.385} & \text{if } H_{\text{w, left 2D}} > H_{\text{w, left 1D}} > 0 \end{cases} \quad (30) $$

where $q_{\text{in, left}}$ is the discharge evaluated on the left side of the rill, $H_{\text{w, left 2D}} = H_{\text{left 2D}} - H_{\text{Crest, left}}$ and $H_{\text{w, left 1D}} = H_{\text{1D}} - H_{\text{Crest, left}}$, with $H_{\text{Crest, left}}$ being the ridge height of the rill on the left, and $C_d$ a user-selected discharge coefficient, the value of which depends on the shape of the weir (Henderson, 1966). This value is the same for all cases in the present work, being $C_d = 0.6$ (Echeverribar et al., 2019). Figs. 5 and 6 show the most important variables of a generic case. In the same way as in (30), the unit flow rate out of the rill is:

$$ q_{\text{out, left}} = \begin{cases} 0 & \text{if } H_{\text{w, left 2D}} \le 0, H_{\text{w, left 1D}} \le 0 \\ \frac{2}{3} \sqrt{2g} C_d H_{\text{w, left 1D}}^{3/2} & \text{if } H_{\text{w, left 1D}} > 0, H_{\text{w, left 2D}} \le 0 \\ \frac{2}{3} \sqrt{2g} C_d H_{\text{w, left 1D}}^{3/2} \left(1.0 - \left(\frac{H_{\text{w, left 2D}}}{H_{\text{w, left 1D}}}\right)^{3/2}\right)^{0.385} & \text{if } H_{\text{w, left 1D}} > H_{\text{w, left 2D}} > 0 \end{cases} \quad (31) $$

Analogous to the Eqs. (30) and (31), the discharges $q_{\text{in, right}}$ and $q_{\text{out, right}}$ for the right side of the rill are obtained. Thus, the discharges $q_{\text{in}}$ and $q_{\text{out}}$ of the expressions (7) and (9) are $q_{\text{in}} = q_{\text{in, left}} + q_{\text{in, right}}$ and $q_{\text{out}} = q_{\text{out, left}} + q_{\text{out, right}}$.

For reasons of numerical stability, the volume transfer depends on the values reached by the water depth with respect to the critical depth (Echeverribar et al., 2019):

* In most cases, when the water depth at both cells (1D and 2D) exceeds the critical water height $h_c = (q^2g^{-1})^{1/3}$, a transfer volume on the left side of the rill is established for a 1D cell i as:

$$ V_{\text{transf } i, \text{ left 1D}} = (q_{\text{in, left}} - q_{\text{out, left}}) B \Delta t_{2D} \Delta s_i \quad (32) $$

where $\Delta t_{2D}$ is the time step of the 2D numerical scheme at the instant of transfer and $\Delta s_i$ is the i cell size of the 1D model, here assumed to be equal to the length of the 2D cell. The value of $V_{\text{transf } j, \text{ left 2D}}$ is the same as in the 1D cell but with opposite sign.

* If the water depth $h$ at one of the 2D cells or the 1D rill cell exceeds the critical water depth $h_c$ associated with the discharge

<figure><img src="image_11.png" /><figcaption>Fig. 8. Schematic of the numerical calculation using the rill boundary condition.</figcaption></figure>

$q = |q_{\text{in}} - q_{\text{out}}|$ being transferred, the volume transport method (VTM, Echeverribar et al., 2019) is used. The critical water depth $h_c$ is imposed in the cell i where $h_i < h_c$, assuming an added volume in that cell. To ensure mass conservation, the following volume must be subtracted in the cell j with the highest water depth:

$$ V_{\text{transf}} = (h_c - h_i) S_i \quad (33) $$

where $S_i$ is the area of the cell i ($S_i = \Omega_i$ for a 2D cell and $S_i = B \Delta s_i$ for a 1D cell).

The difference between the case with volume entering the rill and the case with volume leaving it is, fundamentally, the sign of the transfer.

### 4.2.1. Mass conservation

In case of uniform depth conditions it is easy to verify that our approach is fully conservative. The fundamental aspect that allows coupling the 2D model with the 1D model through the use of a 1D rill internal boundary condition is the mass conservation at the moment of volume transfer from one model to the other. Total mass must be accounted for in both models in order to avoid mass conservation errors. Hence, at the moment of transfer from the 2D model to the 1D model or vice versa, certain situations may occur in which mass is not well accounted for in the model, which may produce a mass
<table><caption>**Table 1**
Computational results for test cases.</caption><thead><tr><th>Test</th><th>N<sub>cells</sub> 2D</th><th>Computational cost 2D (s)</th><th>N<sub>cells</sub> 2D-1D</th><th>Computational cost 2D-1D (s)</th></tr></thead><tbody><tr><td>1</td><td>4023</td><td>500.38</td><td>201</td><td>1.33</td></tr><tr><td>2</td><td>1770</td><td>60.40</td><td>181</td><td>0.21</td></tr><tr><td>3</td><td>110 912</td><td>18 140</td><td>4420</td><td>7.90</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Test case 1.</figcaption></figure>

conservation error. Therefore, the input and output volume of the models must be accurately accounted for. Thus, to conserve mass in both models, the volumes in (32) and (33) will be used:

$$V_{1D}^{n+1} = \sum_{i=1}^{N_{1D}} h_i^{n+1} B \Delta s_i = V_{1D}^{n} \pm V_{\text{transf}} \quad (34)$$

where $V_{1D}^n$ is the volume in the 1D model before the volume transfer between the models and $N_{1D}$ is the number of 1D cells. $V_{\text{transf}} = V_{\text{transf, left}} + V_{\text{transf, right}}$ depends on the case and its sign is determined by whether the volume is inlet (–) or outlet (+) to the rill. For the 2D model, it will be set as total volume:

$$V_{2D}^{n+1} = \sum_{i=1}^{N_{2D}} h_i^{n+1} \Omega_i = V_{2D}^{n} \mp V_{\text{transf}} \quad (35)$$

where $V_{\text{old 2D}}$ is the volume held in the 2D model before volume transfer between models and $N_{2D}$ is the number of 2D cells. The sign of $V_{\text{transf}}$ must be opposite of that Eq. (34) to ensure mass conservation.

### 4.3. Network of several rills (junction)

One of the possibilities of the internal boundary condition is to represent networks of rills, formed by gullies in a watershed, by several tributary streams (Burguete et al., 2014). Therefore, a junction strategy must be developed in order to obtain a smooth evolution of the conserved variables at the junction region.

The volume variation in the junction region during a time step is:

$$V_{J}^{n+1} - V_{J}^{n} = \left( \sum_{i=1}^{N_{\text{in}}} q_{i}^{n} - \sum_{i=1}^{N_{\text{out}}} q_{i}^{n} \right) B \Delta t_{2D} \quad (36)$$

where $V_J$ is the total volume of the rills junction, $N_{\text{in}}$ is the number of cells that introduce discharge in the junction (inlet cells) and $N_{\text{out}}$ is the number of cells that draw flow from the junction (outlet cells). Finally, writing the current volume at the junction $V_J^n = \sum_{i=1}^{N_{\text{in}}+N_{\text{out}}} h_i^n B \Delta s_i$:

$$V_{J}^{n+1} = \sum_{i=1}^{N_{\text{in}}+N_{\text{out}}} h_{i}^{n} B \Delta s_{i} + \left( \sum_{i=1}^{N_{\text{in}}} q_{i}^{n} - \sum_{i=1}^{N_{\text{out}}} q_{i}^{n} \right) B \Delta t_{2D} \quad (37)$$

Using (37), a new level ($H_J^{n+1}$) is obtained for the next time step. Thus, a distinction is made between inlet and outlet cells:

*   **Inlet cell:** The obtained level $H_J^{n+1}$ is imposed as long as its value is higher than the value of the cell.

*   **Outlet cell:**

– If, for an outlet cell $i$, $q_i \neq 0$ and $h_i > 10^{-6}$ m, the discharge:

$$q_i^{n+1} = q_{\text{total in}} = \frac{\sum_{i=1}^{N_{\text{in}}} q_i^n}{N_{\text{out, 0}}} \quad (38)$$

is imposed in the cell $i$ for the next time $n+1$, where $N_{\text{out, 0}}^n$ is the number of outlet cells with non-zero discharge at time $n$. For stability reasons, if $h_i \le 10^{-6}$ m, the next conditions are imposed:

$$h_i^{n+1} = q_{\text{total in}} \Delta t_{2D}^n, \quad q_i^{n+1} = 0 \quad (39)$$

In the case of having a supercritical flow (Fig. 7), a further condition must be imposed in addition to the discharge in the case of having a water height above the tolerance. This condition is the one in the inlet cell, i.e., the level $H_J^{n+1}$ is also imposed.

– If, on the other hand, $q_i = 0$, the level $H_J^{n+1}$ is imposed on the cell $i$.

### 4.4. Numerical calculation with both models

The union between the one-dimensional numerical scheme and the two-dimensional numerical scheme is performed following the scheme in Fig. 8. In it, it can be seen that first the two-dimensional model is calculated, obtaining a time interval $\Delta t_{2D}^n$ at time step $n$. After the two-dimensional calculation, the one-dimensional model is calculated, and therefore, the hydrodynamics in the grooves, obtaining a time interval $\Delta t_{1D}^n$. For both models to advance synchronously, we have the following situations:

*   If $\Delta t_{1D}^n > \Delta t_{2D}^n$, $\Delta t_{1D}^n = \Delta t_{2D}^n$ is imposed and the one-dimensional hydrodynamics is calculated with this step, since there are no numerical stability problems when using a time step smaller than the one obtained.

*   If $\Delta t_{1D}^n = \Delta t_{2D}^n$, one-dimensional hydrodynamics is calculated with this step.

*   If $\Delta t_{1D}^n < \Delta t_{2D}^n$, one-dimensional hydrodynamics is calculated until $\sum_i \Delta t_{1D, i}^n = \Delta t_{2D}^n$ so that both models advance synchronously, always choosing the time interval $\Delta t_{1D, i} = \min(\Delta t_{1D, i}, \Delta t_{2D}^n - \sum_i \Delta t_{1D, i}^n)$. This implies a variation of the 1D time step for the final iteration, that will be computed as the minimum between the 1D time step and the remaining time between the 2D time step and the sum of the 1D time step for all the iterations performed. However, as the literature says (Murillo and García-Navarro, 2010), the imposition of a time step smaller than the calculated one will not produce numerical instabilities.

## 5. Validation

A series of test cases are presented to validate the correct operation of the rill internal boundary condition. In the test cases, neither precipitation nor infiltration is considered, so $f = R = 0$.

### 5.1. Test case 1

A 10 m × 10 m square domain with no slope (Fig. 9) is first assumed. It contains a rill of length $L = 8$ m aligned with the x axis of rectangular
<figure><img src="image_13.png" /><figcaption>Fig. 10. 2D mesh (a) and 2D mesh with rill (b) for the test case 1.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Temporal evolution of the water surface level (a) and the velocity module (b) at the probe location for the test case 1.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 12. Longitudinal profile along the x coordinate at the y coordinate y = 5.0 m for 3 different states.</figcaption></figure>

cross section formed by a width $B = 0.2$ m and a depth $z_{min} = 0.2$ m. Flow is assumed in the y direction. The initial conditions are dry bed everywhere. The inlet hydrograph imposed as the upstream boundary condition goes from zero and grows steadily to $0.015 \text{ m}^3\text{s}^{-1}$ at the end of the simulation. A uniform flow boundary condition is imposed at the downstream section. The Manning coefficient is $n = 0.032 \text{ s m}^{-1/3}$. A probe is located between the end of the rill at $P = (5.00; 4.50)$ m, in order to compare the results provided by the fully 2D mesh (Fig. 10a) and the mesh with the rill (Fig. 10b), which is modeled by 1D cells with $\Delta s = 0.73$ m. The results of the temporal evolution of the water surface level and the velocity module are shown in Fig. 11. Moreover, the computational cost for both meshes is shown in Table 1. As it can be seen in Fig. 11, the meshes provide almost identical results in accuracy, but their computational times are really different as it can be seen in Table 1, being the 2D-1D mesh more computational efficient. In order to see the influence that the rill has on the spatial evolution of the flow, a longitudinal profile is added along the x-coordinate at $y = 5.0$ m. This profile is shown in Fig. 12, where the spatial evolution of the water surface level of the 2D model and the 2D-1D model along the profile can be observed for 3 different states: State 0 ($t = 0.0$ s) where the rill is completely dry, State 1 ($t = 360.0$ s) where the rill has just filled, and State 2 ($t = 540.0$ s) where considerable time has passed since the rill has filled. As it can be seen in Fig. 12, the differences between the two models are primarily until the rill is filled. Thereafter, the results obtained by both models are practically the same, noting that the discretization of the rill using 2D cells or its modeling by means of an internal boundary condition does not entail appreciable differences. To compare the results between both meshes, the root-mean-square error (RMSE) (Barnston, 1992):

$$ \text{RMSE} = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (x_i - \hat{x}_i)^2} \qquad (40) $$

is used, where $x_i$ are the results at the probe location provided by the 2D-1D mesh for a dumping time $i$, $\hat{x}_i$ are the results at the probe location provided by the 2D mesh for a dumping time $i$, and $N$ is the total number of dumping times. For this case, the root-mean-square error is $\text{RMSE} = 1.39 \cdot 10^{-4}$ m for the water surface level and $\text{RMSE} = 8.77 \cdot 10^{-4}$ m/s for the velocity module.

### 5.1.1. Analysis of the rill width

In order to verify the influence that the size of the rill has, the width of the proposed rill is modified to verify how this aspect affects the results obtained. In this way, the same previous case is simulated with two different widths: a smaller width ($B = 0.1$ m) and a larger width ($B = 0.4$ m). The initial condition, the boundary conditions and the rest of the rill characteristics are the same as in the previous case. The results with $B = 0.1$ m are shown in Fig. 13 and the results with $B = 0.4$ m are shown in Fig. 14. The root-mean-square error for the test case with $B = 0.1$ m is RMSE = $9.55 \cdot 10^{-5}$ m for the water surface level and RMSE = $5.50 \cdot 10^{-4}$ m/s for the velocity module. For the test case with $B = 0.4$ m, the root-mean-square error is RMSE = $1.14 \cdot 10^{-3}$ m for the water surface level and RMSE = $1.15 \cdot 10^{-2}$ m/s for the velocity module. As it can be seen in Figs. 13 and 14, and comparing the values of RMSE, the increase in the width of the rills implies a lower accuracy in the results because the size of the rill starts to be non-negligible compared to the dimensions of the domain.

### 5.2. Test case 2

The second test case is the same dry 10 m × 10 m square domain with no slope (Fig. 15). In this case, a rill of length $L = 3$ m, width $B = 0.2$ m and depth $z_{min} = 0.2$ m is assumed aligned with the y axis. Flow is imposed along the y axis according to a hydrograph from zero to $0.015 \text{ m}^3\text{s}^{-1}$ at the upstream boundary condition and a uniform flow is imposed as the downstream boundary condition. The Manning coefficient is $n = 0.032 \text{ s m}^{-1/3}$. A probe is again located between the end of the orifice and the downstream boundary condition, at $P = (5.00; 3.10)$ m, in order to compare the results provided by the fully 2D mesh (Fig. 16a) and the mesh with the rill as IBC (Fig. 16b), which is modeled by 1D cells with $\Delta s = 0.6$ m. The temporal evolution of the water surface level and the velocity module are shown for both meshes in Fig. 17 and the computational results in Table 1. The results of both meshes are really similar for the water surface level, while for the velocity module the difference between the models is more appreciable. Moreover, the computational cost for 2D mesh is higher than for 2D-1D mesh as it can be seen in Table 1. Using Eq. (40), $\text{RMSE} = 6.11 \cdot 10^{-4}$ m for water surface level and $\text{RMSE} = 5.09 \cdot 10^{-3}$ m/s for the velocity module are obtained.
<figure><img src="image_16.png" /><figcaption>Fig. 13. Temporal evolution of the water surface level (a) and the velocity module (b) at the probe location for the test case 1 with $B = 0.1$ m.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. Temporal evolution of the water surface level (a) and the velocity module (b) at the probe location for the test case 1 with $B = 0.4$ m.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Test case 2.</figcaption></figure>

### 5.3. Test case 3

The third test case is a 100 m × 10 m 2D domain with no slope (Fig. 18). In this case, it contains a T-shaped orifice with a major side of 40 meters aligned with the x axis, a minor side of 7 m, a width $B = 0.2$ m and a depth $z_{min} = 0.2$ m. Flow is assumed along the x axis starting again from dry conditions and corresponding to a hydrograph from zero to $0.05 \ m^3s^{-1}$ imposed as the upstream boundary condition

and a uniform flow is imposed as the downstream boundary condition. The Manning coefficient is $n = 0.032 \ s \ m^{-1/3}$. The IBC involves 3 rills and a junction. A probe is located between the end of the orifice and the downstream boundary condition at $P = (60.50; 5.00)$ m, in order to compare the results provided by the fully 2D mesh (Fig. 19) and the mesh with 3 rills (Fig. 20), which are modeled by 1D cells with $\Delta s = 0.73$ m. The results provided by both meshes are shown in Fig. 21. Only small differences when the rill is full and the flow starts to come out of it can be observed for the water depth (see Fig. 21a). However, the difference between 2D and 2D-1D models in the velocity module is higher, as it can be seen in Fig. 21b. Table 1 shows that the computational cost of the 2D mesh is much higher than that of the 2D mesh with IBC. Using (40), RMSE $= 2.17 \cdot 10^{-3}$ m for the water surface level and RMSE $= 8.73 \cdot 10^{-3}$ m/s for the velocity module are obtained for this test case.

## 6. Application to the Arnás river catchment

The Arnás river catchment is located in the Borau valley (Northern Spanish Pyrenees). In the last decades, an exhaustive experimental study of the water-table dynamics, sediment transport and erosion processes, runoff generation and streamflow response has been carried out (Lana-Renault et al., 2007; García-Ruiz et al., 2005). From a numerical perspective, the physical processes of rainfall, infiltration, runoff
<figure><img src="image_19.png" /><figcaption>Fig. 16. 2D mesh (a) and 2D mesh with rill (b) for the test case 2.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 17. Temporal evolution of the water surface level (a) and the velocity module (b) at the probe location for the test case 2.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 18. Test case 3.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 19. 2D mesh for the test case 3.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 20. 2D mesh with 3 rills for the test case 3.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 21. Temporal evolution of the water surface level (a) and the velocity module (b) at the probe location for the test case 3.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 22. 3D representation of the Arnás catchment hypsometry map.</figcaption></figure>
<table><caption>**Table 2**
Number of cells for all the assessed meshes.</caption><thead><tr><th></th><th>Mesh 0</th><th>Mesh 1</th><th>Mesh 2</th><th>Mesh 3</th><th>Mesh 4</th><th>Mesh 5</th></tr></thead><tbody><tr><th>N<sub>cells</sub></th><td>404 152</td><td>166 573</td><td>113 874</td><td>92 311</td><td>61 940</td><td>38 553</td></tr></tbody></table>

generation, groundwater flow, erosion and sedimentation have been simulated in detail using different methods (López-Barrera et al., 2011; Fernández-Pato et al., 2016). The main goal of this case (see Fig. 22) is to assess if the rill internal boundary condition can outperform a fully 2D mesh in a realistic application of rainfall runoff flow at a mountain basin level where several previous analysis have been performed (López-Barrera et al., 2011; Fernández-Pato et al., 2016; Caviedes-Voullième et al., 2012). The first step is to find the coarsest mesh that allows the full 2D model to generate accurate numerical results. Starting from a simulation carried out with an ultra-fine reference mesh (Mesh 0) of rainfall over dry terrain, the number of cells has been progressively reduced and the hydrographs at the outlet of the basin have been compared. The refinement of the reference mesh (Mesh 0) has been carried out taking into account the maximum refinement level imposed by the resolution of the available digital terrain model of the basin (2 m × 2 m). Table 2 shows the number of cells ($N_{cells}$) for all the assessed meshes.

A single rainfall event with a duration of 2.4 days is considered for this test. Fig. 23 shows the measured hyetograph, registered by a rain gauge with a frequency of 1 h. All the domain boundaries are closed except for the outlet point, located at the position of the flowmeter (see Fig. 22). Regarding the Manning roughness coefficient, a uniform value of $n = 0.05 \text{ s m}^{-1/3}$ is set for the sake of simplicity. The infiltration parameters are $k = 7 \cdot 10^{-4} \text{ s}^{-1}$, $f_c = 2 \cdot 10^{-6} \text{ m/s}$ and $f_0 = 2.2 \cdot 10^{-5} \text{ m/s}$, following the values used in previous works (Fernández-Pato et al., 2016).

The outlet numerical hydrographs for all the assessed meshes are compared in Fig. 24. Several conclusions are reached here. On the one hand, all the meshes except Mesh 5 are capable of capturing the main peak of the outlet hydrograph, so this mesh is discarded. The hydrographs corresponding to Mesh 0 and Mesh 1 are practically identical so it is clear that the coarsening of Mesh 1 is not enough. The two secondary peaks of the final part of the hydrograph are correctly simulated with any of the meshes considered (see Fig. 26). On the other hand, paying attention to the first stages of the event, we see significant differences among meshes (see Fig. 25). Differences up to 50% in the third discharge peak automatically discard Mesh 4, while Mesh 3 fails at reproducing adequately the time at which the hydrograph begins. Taking all these considerations into account, we can conclude that the optimal mesh for the 2D model in this particular case is Mesh 2 with 113874 computational cells.

When the optimum 2D mesh is obtained, the different regions where the main gullies are located in the domain will be modeled using rills as IBC. These regions must be discretized correctly because of their complexity and their irregularity (Fig. 22), therefore, the size of the cells in this zones must be smaller than the rest of the domain (Fig. 27a). As shown in Fig. 28, there are 7 main gullies, but two of them are much smaller than the rest. For this reason, two meshes are proposed, one with the 5 bigger rills and another with 7 rills. In this particular case, taking into account the topography of the domain, the secondary rills (2 to 7) are not connected to the main one because the mouth of the secondary rills in the main rill is an extension of the width of the gullies where the water is accumulated, so modeling this region using a 1D rill would be a high error. The main difference between the fully 2D mesh (M2D) and the 2D mesh with rills (M2D5R and M2D7R) is the size of cells in the gullies, as it can be seen in Fig. 27 and in Table 3. To compare the different meshes, a probe is placed at the position where the flowmeter and turbidimeter are located, very close

<figure><img src="image_26.png" /><figcaption>**Fig. 23.** Hyetograph for the Arnás catchment test case.</figcaption></figure>

<table><caption>**Table 3**
Computational results for Arnás case.</caption><thead><tr><th></th><th>2D Mesh</th><th>2D Mesh + 5 rills</th><th>2D Mesh + 7 rills</th></tr></thead><tbody><tr><th>N<sub>cells</sub></th><td>113 874</td><td>78 038</td><td>77 950</td></tr><tr><th>Computational cost (s)</th><td>32 385.76</td><td>8970.98</td><td>8924.24</td></tr><tr><th>Speed-Up</th><td>-</td><td>3.61</td><td>3.63</td></tr></tbody></table>

<table><caption>**Table 4**
RMSE for Arnás case.</caption><thead><tr><th></th><th>2D Mesh + 5 rills</th><th>2D Mesh + 7 rills</th></tr></thead><tbody><tr><th>RMSE (without correction)</th><td>0.170 m<sup>3</sup>/s</td><td>0.176 m<sup>3</sup>/s</td></tr><tr><th>RMSE (with correction)</th><td>0.158 m<sup>3</sup>/s</td><td>0.172 m<sup>3</sup>/s</td></tr></tbody></table>

to the downstream boundary condition. The temporal evolution of the outflow discharge is shown in Fig. 29 and the computational cost for each mesh is shown in Table 3. The three meshes provide very similar results for most of the event, except for the first rain period, where M2D5R and M2D7R provide a significantly higher outflow discharge than the M2D. This may be due to the fact that, as the cell size increases near the gullies, the infiltration rate in these regions is lower than for M2D, which initially results in a lower infiltration rate, which over time balances out in all 3 meshes. For this reason, the initial infiltration rate $f_0$ is slightly increased to $f_0 = 2.56 \cdot 10^{-5} \text{ m/s}$ in the meshes with rills, obtaining the temporal evolution of the outflow discharge shown in Fig. 30. This figure shows a greater similarity between the results obtained with each mesh. Figs. 31, 32 and 33 are added to show the water depth in the entire domain using both models, showing great similarity in the results shown. Moreover, the computational costs are shown in the Table 3, where it can be seen that the meshes with rills are much more computationally efficient than the 2D mesh by significantly not only reducing the number of cells, but also increasing the time step size that decreases when the grid size gets smaller (leading to a larger number of time steps to complete the simulation). The RMSE (Eq. (40)) is obtained for the cases with and without infiltration correction. The values of RMSE are shown in Table 4. From these results it can be concluded that the use of 5 rills provides information more similar to that provided by the 2D mesh. In addition, the infiltration correction reduces the error in the meshes with rills, being the mesh with 5 rills and infiltration correction the mesh with the lowest error.
<figure><img src="image_27.png" /><figcaption>Fig. 24. Outlet hydrographs for all meshes.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 25. Zoom view of the beginning of the outlet hydrographs.</figcaption></figure>

## 7. Conclusions

In this paper, a rill internal boundary condition for narrow regions on a 2D mesh is presented, in order to reduce the computational cost of realistic event simulations without losing the accuracy offered by fully 2D meshes.

The main difference with previous works is the simplicity of this coupling, where the one-dimensional cells are not physical cells in the mesh, but their presence is shown in the mesh with an increase of 2D cells in these regions, with loss of topographic information that is compensated by the presence of the rill boundary condition in these regions.

Three validation cases are presented to test the accuracy and computational cost offered by this internal boundary condition. In all three cases, very similar results for the water surface level are obtained between the use of a fully 2D mesh and the use of rills to model the narrow regions of the domain. This accuracy is reflected by the low RMSE values obtained. However, the difference between the models

are higher when the velocity module results are compared. The main reason of this difference is the discharge used in the 1D rill, which models the flow entering the rill as if it were a weir. Therefore, some error is introduced in this approximation, although, as can be seen in the RMSE values obtained for the velocity module, this error is still low. In addition, the computational times offered by meshes with rills are much lower than those offered by 2D meshes, so that the use of this type of boundary condition in test cases optimizes the mesh significantly with practically no loss of accuracy in the results.

Subsequently, a realistic case, the Arnás river catchment, is studied in order to check if this internal boundary condition can be useful for certain regions with gullies or narrow regions that significantly increase a 2D mesh by requiring cells of a reduced size to capture the topographic information of these areas correctly. Due to the presence of numerous gullies in this basin, two meshes with different number of rills are developed, one with 5 rills and another with 7 rills. The results obtained with these two meshes are similar to those offered
<figure><img src="image_29.png" /><figcaption>Fig. 26. Zoom view of the end of the outlet hydrographs.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. 27. 2D mesh without rill (a) and with rill (b) for the Arnás case.</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>Fig. 28. Arnás case domain with the different regions where the rills are located.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Fig. 29. Temporal evolution of the discharge for 2D mesh, 2D mesh with 5 rills and 2D mesh with 7 rills, without the infiltration correction.</figcaption></figure>

by a completely 2D mesh, except in the first hours of raining, where the results differ notably, as it can be seen in Fig. 29. This aspect may be due to the fact that, since there are fewer cells in the areas where the rills are located, these meshes do not initially infiltrate as much water as the 2D mesh does. For this reason, the initial infiltration rate is slightly increased in these two meshes. In this way, a greater similarity in the results is observed (see Fig. 30), reducing the RMSE value in comparison with the meshes without infiltration correction (see Table 4). With respect to computational times, the use of rills represents a great reduction in computational cost as it can be seen in Table 3.

Thus, it can be concluded that the implemented internal boundary condition is an accurate and efficient tool for modeling narrow regions in large domains, which significantly increase the number of cells, and therefore, the computational cost. It is important to note, however, that in a catchment, the approach requires prior knowledge of the location of the rills and therefore a detailed understanding of hydrographic network. Moreover, the rill width should be negligible at the domain scale, otherwise some inaccuracies can be observed (see Figs. 13 and 14). When this information is not available, the approach seems inapplicable unless a preliminary simulation is conducted with a very fine computational mesh to identify the hydrographic network.
<figure><img src="image_33.png" /><figcaption>Fig. 30. Temporal evolution of the discharge for 2D mesh, 2D mesh with 5 rills and 2D mesh with 7 rills, with the infiltration correction.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. 31. Spatial distribution of the water depth for 2D mesh without rill (a) and with rills (b) for the Arnás case in *t* = 2.25 hours.</figcaption></figure>
<figure><img src="image_35.png" /><figcaption>Fig. 32. Spatial distribution of the water depth for 2D mesh without rill (a) and with rills (b) for the Arnás case in $t = 3.25$ hours.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Fig. 33. Spatial distribution of the water depth for 2D mesh without rill (a) and with rills (b) for the Arnás case in $t = 4.25$ hours.</figcaption></figure>

### CRediT authorship contribution statement

P. Vallés: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. J. Fernández-Pato: Writing – original draft, Validation, Software. M. Morales-Hernández: Supervision, Methodology, Investigation, Formal analysis, Conceptualization. I. Echeverría: Software, Methodology. P. García-Navarro: Writing – review & editing, Writing – original draft, Supervision, Project administration, Investigation, Funding acquisition, Formal analysis, Conceptualization.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: All

authors report financial support was provided by Aragón Government, DGA (Spain).

### Data availability

No data was used for the research described in the article.

### Acknowledgments

This work was partially funded by the Aragón Government, DGA (Spain), through the Fondo Europeo de Desarrollo Regional, FEDER. Moreover, the authors would like to thank HYDRONIA LLC for their help during the development of this work with the motivation to enhance the software RiverFlow2D.
## References

Aggett, G., Wilson, J., 2009. Creating and coupling a high-resolution DTM with a 1-D hydraulic model in a GIS for scenario-based assessment of avulsion hazard in a gravel-bed river. Geomorphology 113 (1), 21–34. http://dx.doi.org/10.1016/j.geomorph.2009.06.034. Understanding earth surface processes from remotely sensed digital terrain models.

Arcement, G., Schneider, V., 1984. Guide for Selecting Manning's Roughness Coefficients for Natural Channels and Flood Plains. In: U.S. Geological Survey. Water-supply paper, (2339).

Barnston, A.G., 1992. Correspondence among the correlation, RMSE, and Heidke forecast verification measures; Refinement of the Heidke score. Weather Forecast. 7 (4), 699–709. http://dx.doi.org/10.1175/1520-0434(1992)007<0699>CATCRA>2.0.CO;2.

Bates, P.D., 2000. Development and testing of a subgrid-scale model for moving-boundary hydrodynamic problems in shallow water. Hydrol. Process. 14 (11–12), 2073–2088. http://dx.doi.org/10.1002/1099-1085(20000815/30)14:11/12<2073::AID-HYP55>3.0.CO;2-X.

Bcc Waterways Program Staff, Witheridge, G., Brisbane (Qld.). City Council, 2000. Natural Channel Design Guidelines. Appendix C. In: Technical Document, Brisbane City Council.

Betsholtz, A., Nordlöf, B., 2017. Potentials and limitations of 1D, 2D and coupled 1D-2D flood modelling in HEC-RAS.

Bomers, A., Schielen, R.M.J., Hulscher, S.J.M.H., 2019. The influence of grid shape and grid size on hydraulic river modelling performance. Environ. Fluid Mech. 1273–1294. http://dx.doi.org/10.1007/s10652-019-09670-4.

Briganti, R., Dodd, N., 2009. Shoreline motion in nonlinear shallow water coastal models. Coast. Eng. 56 (5), 495–505.

Brunner, G.W., 2014. Combined 1D and 2D modeling with HEC-RAS. US Army Corps Eng.

Burguete, J., Lacasta, A., García-Navarro, P., 2014. SURCOS: A software tool to simulate irrigation and fertigation in isolated furrows and furrow networks. Comput. Electron. Agric. 103, 91–103. http://dx.doi.org/10.1016/j.compag.2014.02.004.

Casulli, V., 2009. A high-resolution wetting and drying algorithm for free-surface hydrodynamics. Internat. J. Numer. Methods Fluids 60 (4), 391–408. http://dx.doi.org/10.1002/fld.1896.

Caviedes-Voullième, D., García-Navarro, P., Murillo, J., 2012. Influence of water structure on 2D full shallow water equations and SCS Curve Number simulation of rainfall/runoff events. J. Hydrol. 448–449, 39–59. http://dx.doi.org/10.1016/j.jhydrol.2012.04.006.

Costabile, P., Costanzo, C., Macchione, F., 2017. Performances and limitations of the diffusive approximation of the 2-d shallow water equations for flood simulation in urban and rural areas. Appl. Numer. Math. 116, 141–156. http://dx.doi.org/10.1016/j.apnum.2016.07.003. New Trends in Numerical Analysis: Theory, Methods, Algorithms and Applications (NETNA 2015).

Costabile, P., Macchione, F., Matale, L., Petaccia, G., 2015. Flood mapping using LIDAR DEM. Limitations of the 1-D modeling highlighted by the 2-D approach. Nat. Hazards 77 (2), 181–204.

Dazzi, S., Vacondio, R., Mignosa, P., 2020. Internal boundary conditions for a GPU-accelerated 2D shallow water model: Implementation and applications. Adv. Water Resour. 137, 103525. http://dx.doi.org/10.1016/j.advwatres.2020.103525.

Defina, A., 2000. Two-dimensional shallow flow equations for partially dry areas. Water Resour. Res. 36 (11), 3251–3264. http://dx.doi.org/10.1029/2000WR090167.

Dhondia, J., Stelling, G., 2002. Application of one dimensional - two dimensional integrated hydraulic model for flood simulation and damage assessment. In: Falconer, R., Lin, B., Harris, E., Wilson, C. (Eds.), Proceedings of the Fifth International Conference on Hydroinformatics Model Development and Data Management. International Water Association (IWA), United Kingdom, pp. 265–282, Fifth international conference on hydroinformatics, Cardiff, UK ; Conference date: 01-07-2002 Through 05-07-2002.

Echeverribar, I., Morales-Hernández, M., Brufau, P., García-Navarro, P., 2019. Use of internal boundary conditions for levees representation: application to river flood management. Environ. Fluid Mech. 1253–1271. http://dx.doi.org/10.1007/s10652-018-09658-6.

Fernández-Pato, J., Caviedes-Voullième, D., García-Navarro, P., 2016. Rainfall/runoff simulation with 2D full shallow water equations: Sensitivity analysis and calibration of infiltration parameters. J. Hydrol. 536, 496–513. http://dx.doi.org/10.1016/j.jhydrol.2016.03.021.

Finaud-Guyot, P., Delenne, C., Guinot, V., Llovel, C., 2011. 1D-2D coupling for river flow modeling. C. R. Méc. 339 (4), 226–234. http://dx.doi.org/10.1016/j.cme.2011.02.001.

García-Ruiz, J., Arnáez, J., Beguería, S., Seeger, M., Martí-Bono, C., Regüés, D., Lana-Renault, N., White, S., 2005. Runoff generation in an intensively disturbed, abandoned farmland catchment, Central Spanish Pyrenees. CATENA 59 (1), 79–92. http://dx.doi.org/10.1016/j.catena.2004.05.006.

Gejadze, I., Monnier, J., 2007. On a 2D 'zoom' for the 1D shallow water model: Coupling and data assimilation. Comput. Methods Appl. Mech. Engrg. 196 (45), 4628–4643. http://dx.doi.org/10.1016/j.cma.2007.05.026.

Henderson, F., 1966. Open Channel Flow. In: Macmillan Series in Civil Engineering, Macmillan.

Horritt, M.S., Bates, P.D., 2002. Evaluation of 1D and 2D numerical models for predicting river flood inundation. J. Hydrol. 268, 89–99.

Horton, R.E., 1933. The Rôle de infiltration in the hydrologic cycle. EOS Trans. Am. Geophys. Union 14 (1), 446–460. http://dx.doi.org/10.1029/TR014i001p00446.

Hu, R., Fand, F., Salinas, P., Pain, C.C., Sto. Domingo, N.D., Mark, O., 2019. Numerical simulation of floods from multiple sources using an adaptive anisotropic unstructured mesh method. Adv. Water Resour. 123, 173–188.

Kitts, D., Syme, W., Gao, S., Colleutt, G., Ryan, P., 2020. Mesh orientation and cell size sensitivity in 2D SWE solvers. In: RiverFlow 2020. Taylor & Francis Group, London, p. 9.

Lana-Renault, N., Latron, J., Regüés, D., 2007. Streamflow response and water-table dynamics in a sub-Mediterranean research catchment (Central Pyrenees). J. Hydrol. 347 (3), 497–507. http://dx.doi.org/10.1016/j.jhydrol.2007.09.037.

LeVeque, R., 1992. Numerical Methods for Conservation Laws. In: Lectures in Mathematics. ETH Zürich, Birkhäuser Basel.

Li, W., Chen, Q., Mao, J., 2009. Development of 1D and 2D coupled model to simulate urban inundation: An application to Beijing Olympic village. Chin. Sci. Bull. 54 (9), 1613–1621. http://dx.doi.org/10.1007/s11434-009-0208-1.

Li, Q., Schlatter, P., Henningson, D.S., 2011. Comparison of SGS models for passive scalar mixing in turbulent channel flows. In: Kuerten, H., Geurts, B., Armenio, V., Fröhlich, J. (Eds.), Direct and Large-Eddy Simulation VIII. Springer Netherlands, Dordrecht, pp. 131–136.

Lin, B., Wicks, J.M., Falconer, R.A., Adams, K., 2006. Integrating 1D and 2D hydrodynamic models for flood simulation. Proc. Inst. Civ. Eng. - Water Manage. 159 (1), 19–25. http://dx.doi.org/10.1680/wama.2006.159.1.19.

López-Barrera, D., Navarro, P.G., Brufau, P., 2011. Sources of uncertainty in the validation of a coupled hydrological-hydraulic simulation model with sediment transport. Houille Blanche 97 (3), 17–22. http://dx.doi.org/10.1051/lhb/201102107.

Marin, J., Monnier, J., 2009. Superposition of local zoom models and simultaneous calibration for 1D-2D shallow water flows. Math. Comput. Simulation 80 (3), 547–560. http://dx.doi.org/10.1016/j.jmatcom.2009.09.001.

Masoero, A., Claps, P., Asselman, N.E.M., Mosselman, E., Di Baldassarre, G., 2013. Reconstruction and analysis of the Po River inundation of 1951. Hydrol. Process. 27 (9), 1341–1348. http://dx.doi.org/10.1002/hyp.9558.

Masood, M., Takeuchi, K., 2012. Assessment of flood hazard, vulnerability and risk of mid-eastern Dhaka using DEM and 1D hydrodynamic model. Nat. Hazards 61 (2), 757–770. http://dx.doi.org/10.1007/s11069-011-0060-x.

Miglio, E., Perotto, S., Saleri, F., 2005. Model coupling techniques for free-surface flow problems: Part I. Nonlinear Anal. TMA 63 (5), 1885–1896. http://dx.doi.org/10.1016/j.na.2005.03.083. Invited Talks from the Fourth World Congress of Nonlinear Analysts (WCNA 2004).

Morales-Hernández, M., García-Navarro, P., Burguete, J., Brufau, P., 2013. A conservative strategy to couple 1D and 2D models for shallow water flow simulation. Comput. & Fluids 81, 26–44.

Murillo, J., García-Navarro, P., 2010. Weak solutions for partial differential equations with source terms: Application to the shallow water equations. J. Comput. Phys. 229, 4237–4368.

Nihei, Y., Ogata, Y., Yoshimura, R., Ito, T., Kashiwada, J., 2024. Subgrid model of water storage in paddy fields for a grid-based distributed Rainfall-Runoff model and assessment of paddy field dam effects on flood control. Water 16 (2), http://dx.doi.org/10.3390/w16020255. URL https://www.mdpi.com/2073-4441/16/2/255.

Petaccia, G., Natale, L., Savi, F., Velickovic, M., Zech, Y., Soares-Frazão, S., 2012. Flood wave propagation in steep mountain rivers. J. Hydroinform. 15 (1), 120–137.

Roe, P., 1981. Approximate Riemann solvers, parameter vectors, and difference schemes. J. Comput. Phys. 43 (2), 357–372. http://dx.doi.org/10.1016/0021-9991(81)90128-5.

Sanders, B., Schubert, J., Detwiler, R., 2010. Parbrezo: a parallel, unstructured grid, godunov type, shallow water code for high resolution flood inundation modeling at the regional scale. Adv. Water Resour. 33 (12), 1456–1467.

Sehili, A., Lang, G., Lippert, C., 2014. High-resolution subgrid models: background, grid generation, and implementation. Ocean Dyn. 64 (4), 519–535. http://dx.doi.org/10.1007/s10236-014-0693-x.

Sotelo, G., 1977. Hidráulica general: fundamentos. Ed. Limusa, México.

Toro, E.F., 1997. Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer, Berlin.

Vacondio, R., Aureli, F., Mignosa, P., Dal Palú, A., 2016. Simulation of the January 2014 flood on the Secchia River using a fast and high-resolution 2D parallel shallow-water numerical scheme. Nat. Hazards 80 (1), 103–125.

Verwey, A., 2001. Latest developments in floodplain modelling - 1D/2D integration. In: 6th Conference on Hydraulics in Civil Engineering: The State of Hydraulics; Proceedings. Barton, A.C.T., URL https://search.informit.org/doi/10.3316/informit.521256778161310.

Villanueva, I., Wright, N.G., 2006. Linking Riemann and storage cell models for flood prediction. Proc. Inst. Civ. Eng. - Water Manage. 159 (1), 27–33. http://dx.doi.org/10.1680/wama.2006.159.1.27.

Vreugdenhil, C.B., 1994. Numerical Methods for Shallow-Water Flow. Kluwer Academic Publishers, Dordrecht, The Netherlands.