# Incorporating the vertical velocity in a coupled Lagrangian–Eulerian approach for particle transport in shallow flows

P. Vallés, J. Segovia-Burillo, M. Morales-Hernández, Volker Roeber, P. García-Navarro

► To cite this version:

P. Vallés, J. Segovia-Burillo, M. Morales-Hernández, Volker Roeber, P. García-Navarro. Incorporating the vertical velocity in a coupled Lagrangian–Eulerian approach for particle transport in shallow flows. Advances in Water Resources, 2025, 205, pp.105085. 10.1016/j.advwatres.2025.105085. hal-05247461

HAL Id: hal-05247461

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

Submitted on 9 Sep 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>

# Incorporating the vertical velocity in a coupled Lagrangian–Eulerian approach for particle transport in shallow flows

P. Vallés<sup>a,b,c</sup>,*, J. Segovia-Burillo<sup>a,c</sup>, M. Morales-Hernández<sup>a,c</sup>, V. Roeber<sup>b,c</sup>,
P. García-Navarro<sup>a,c</sup>

## ARTICLE INFO

Keywords:

Vertical velocity

Shallow water

Finite volume scheme

Debris transport

Lagrangian particle-tracking method

## ABSTRACT

This work presents a method to incorporate vertical velocity into a two-dimensional depth-averaged Shallow Water Equation (2DH SWE) model, thereby improving the calculation of particle trajectories in a Lagrangian Particle Tracking (LPT) framework. The resulting formulation couples Eulerian and Lagrangian approaches. The vertical velocity is also used to modify the dispersion terms in the LPT model. The proposed approximation is first validated—without particle transport—by comparison with Hyperbolic–Elliptic and Hyperbolic-Relaxed Non-Hydrostatic Pressure (NHP) models. The differences between models are minor, confirming the suitability of the vertical velocity approximation for shallow flow problems. Subsequently, the method is applied to particle transport scenarios, demonstrating that including vertical velocity yields more realistic particle trajectories in complex flow situations.

## 1. Introduction

In recent decades, predictive tools based on numerical simulations have been developed and improved, becoming essential for forecasting the evolution of physical systems (Sun et al., 2021; Burton et al., 2020; Ooyama, 1969; Gairola and Bitsuamlak, 2019). Consequently, governments and public institutions have increasingly used these tools to predict and mitigate the negative effects of various natural disasters, such as volcanic eruptions (Poland and Anderson, 2020; Chevrel et al., 2018) and flood events (Thielen et al., 2009; Knijff et al., 2010; Chen et al., 2009). For the latter, the two-dimensional depth-averaged Shallow Water Equations (2DH SWE) models can accurately simulate such disasters with high accuracy and computational efficiency (Sanders et al., 2010; Vacondio et al., 2016). This efficiency arises from the assumption of a hydrostatic pressure distribution and the neglect of vertical velocity. Moreover, implementing these models to enable simulation on GPUs (Castro et al., 2011; Petaccia et al., 2016) or multiple GPUs (Sætra and Brodtkorb, 2012; Xia et al., 2019) has enhanced their efficiency, making them valuable tools for real-time flood forecasting.

The 2DH SWE models are used not only to simulate the evolution of water flow but also to model several geophysical flows, such as mud (Armanini et al., 2009), oil (Echeverribar et al., 2023b), and lava (Ortega-Moya et al., 2024). Furthermore, components like sediments or pollutants can be incorporated into these models using

the flow properties. These elements can be described using either an Eulerian or Lagrangian approach, depending on the problem characteristics (Nordam et al., 2023). In the Eulerian description, sediments or pollutants are modeled as a concentration and transported using expressions based on the advection–diffusion equation (Morgan et al., 1998; Martínez-Aranda et al., 2020). This approach offers a good balance between accuracy and computational efficiency. However, when the number of pollutants is low or when specific processes for each element must be modeled, the Eulerian description may not provide sufficiently detailed results. For this reason, the Lagrangian approach is preferred in such cases (García-Martínez and Flores-Tovar, 1999; Baharvand et al., 2023). The Lagrangian description allows for the modeling of specific processes of the transported elements, such as degradation, fragmentation, and deposition. While this yields more realistic particle trajectories, it also increases computational cost.

Several studies have developed coupled Lagrangian–Eulerian 2DH SWE models for material transport, providing valuable results across many fields (Finaud-Guyot et al., 2023; Persi et al., 2018). However, a key challenge for these coupled models arises when transported particles have a variable vertical position during the simulation. When obstacles are present on the bottom, the depth-averaged horizontal velocities may lead to inaccurate or even nonphysical results. The actual vertical flow velocity component prevents such unrealistic trajectories

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

<footnote><sup>b</sup> E2S–UPPA, SIAME, Université de Pau et des Pays de l'Adour, Anglet, France</footnote>

<footnote>* Corresponding author at: I3A, University of Zaragoza, Zaragoza, Spain.
E-mail address: pvalles@unizar.es (P. Vallés).</footnote>
by either lifting or depositing particles whose vertical positions are lower than the obstacle height. To simulate these effects, 3D models are a possible way to obtain the temporal and spatial evolution of vertical velocity, providing accurate results for Lagrangian particle transport (Jalón-Rojas et al., 2019; Pilechi et al., 2022). However, the involved computational cost reduces the balance between accuracy and efficiency, especially when simulating large-scale problems (Teng et al., 2017).

This work presents an approximation for vertical velocity in a 2DH SWE model to enhance particle trajectories in a Lagrangian transport model developed within the SERGHEI framework (https://gitlab.com/serghei-model/serghei), a multi-dimensional, multi-domain, and multi-physics model for environmental and landscape simulation, designed with an outlook towards Earth system modeling (Caviedes-Voullième et al., 2023). Although similar approximations have been derived in previous works (Ferrari and Saleri, 2004), they have seen limited practical application due to their lack of influence on the evolution of hydraulic variables. In this study, their value is demonstrated in the context of Lagrangian particle transport, where vertical motion might play an important role. To enhance the representation of this motion, we introduce a vertical linear profile for the velocity component, allowing for a more realistic description of particle trajectories while maintaining compatibility with the depth-averaged hydrodynamic model. The vertical velocity is derived from the conservative variables solved by the 2DH SWE model, providing a physically consistent and computationally efficient alternative to fully 3D implementations. Non-Hydrostatic Pressure (NHP) SWE models have been recently developed and refined, providing accurate results in scenarios where vertical velocity plays a significant role (Escalante et al., 2018, 2019; Echeverribar et al., 2023a). Before improving the Lagrangian model with the proposed vertical velocity approximation, its accuracy must be validated by comparing its results with those from two NHP SWE models: a Hyperbolic-Elliptic (HE)-NHP model and a Hyperbolic-Relaxed (HR)-NHP model (Echeverribar et al., 2023a). For simplicity, these comparisons are performed using one-dimensional depth-averaged (1DH) models. Several test cases are simulated to assess the performance of the different models and to determine whether the proposed approach yields comparable vertical velocity values. Once the accuracy of the proposed vertical velocity estimation is evaluated, it is ready to be incorporated into the 2DH SWE model, with the aim to provide more realistic particles trajectories without a relevant increase in computational cost.

The paper is organized as follows: first, the equations used in the 1DH SWE model are presented. Next, for the sake of simplicity, the process for obtaining the vertical velocity approximation is presented in the framework of the 1D system. Subsequently, the numerical scheme is introduced, with a focus on the implementation of the vertical velocity. The equations of the 1D HE-NHP and HR-NHP models are described before the presentation of the test cases. Three analytical test cases are simulated to assess the accuracy of the vertical velocity approximation by comparing with the vertical velocity obtained from the 1D NHP models. Following this, the equations of the 2DH SWE model supplied with the approach to estimate the vertical velocity component and the LPT model are presented, with an emphasis on the vertical velocity approximation and its influence on particle trajectories. The numerical schemes for both the 2DH SWE model and the LPT model are then detailed. Another three analytical test cases are simulated to evaluate the impact of vertical velocity on particle trajectories. Finally, a discussion of the results and the conclusions are provided.

## 2. Governing equations for 1DH model

The 1D dynamic flow equations are presented first with the approach to obtain the vertical velocity. Note that these equations can be derived as a simplification of the 2D Shallow Water Equations (SWEs), under the assumption that the flow has no velocity component in the

transversal (y) direction. The 2DH SWEs themselves are obtained by integrating the incompressible Navier-Stokes equations over the vertical coordinate, assuming hydrostatic pressure distribution and negligible vertical accelerations. This reduction to one dimension is valid when variations in the transversal direction are negligible and the flow is predominantly aligned along the x-axis.

### 2.1. 1DH SWE model

The one-dimensional hydrodynamic model describes the depth-averaged surface flow per unit width by means of the hyperbolic 1DH SWE system based on mass and momentum conservation under hydrostatic pressure hypothesis (Cunge et al., 1989):

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

where:

$$ U = \begin{pmatrix} h \\ h\bar{u} \end{pmatrix}, \quad F = \begin{pmatrix} h\bar{u} \\ (h\bar{u})^2/h + \frac{1}{2}gh^2 \end{pmatrix}, \quad S = \begin{pmatrix} 0 \\ gh(S_0 - S_f) \end{pmatrix} \qquad (2) $$

being h the water depth [L], $\bar{u}$ the x-coordinate depth-averaged velocity flow [$LT^{-1}$], and g the gravitational acceleration [$LT^{-2}$]. $S_0$ represents the bottom slope:

$$ S_0 = -\frac{\partial z_b}{\partial x} \qquad (3) $$

where $z_b$ is the bottom level [L], and $S_f$ represents the friction slope:

$$ S_f = \frac{n^2 |\bar{u}| \bar{u}}{h^{4/3}} \qquad (4) $$

where n is the Manning's roughness coefficient [$TL^{-1/3}$] (Arcement and Schneider, 1984).

### 2.2. Adding the vertical velocity in the 1DH SWE model

For the 1DH SWE model, the system of Eqs. (1) consists of the depth-integrated mass conservation equation and the momentum equation in the x-direction. According to Bristeau et al. (2015), Castro-Orgaz and Hager (2017), an additional equation can be used to approximate the vertical velocity. This equation is derived from the divergence-free condition ($\nabla \cdot \mathbf{v} = 0$, where v denotes the flow velocity vector) and is commonly referred to as the Depth Integrated Incompressibility Condition (DIIC) (Escalante et al., 2018; Yamazaki et al., 2009). The DIIC is obtained by integrating the divergence-free condition over the vertical coordinate from the bottom, $z_b$, to the free surface elevation, $H = h + z_b$:

$$ \int_{z_b}^{H} \frac{\partial u}{\partial x} dz + \int_{z_b}^{H} \frac{\partial w}{\partial z} dz = 0 \qquad (5) $$

Applying the Leibniz rule, the following expressions are obtained for each term from (5):

* For the x-term in (5), the following expression is obtained applying the Leibniz rule:

$$ \int_{z_b}^{H} \frac{\partial u}{\partial x} dz = \frac{\partial}{\partial x} \int_{z_b}^{H} u \, dz - [u]_{H} \frac{\partial H}{\partial x} + [u]_{z_b} \frac{\partial z_b}{\partial x} \qquad (6) $$

and, developing the term $\partial H/\partial x$ considering that $H = h + z_b$:

$$ \int_{z_b}^{H} \frac{\partial u}{\partial x} dz = \frac{\partial}{\partial x} \int_{z_b}^{H} u \, dz - [u]_{H} \left( \frac{\partial h}{\partial x} + \frac{\partial z_b}{\partial x} \right) + [u]_{z_b} \frac{\partial z_b}{\partial x} \qquad (7) $$

Considering the definition of mean horizontal velocity $\bar{u}$:

$$ \bar{u} = \frac{1}{h} \int_{z_b}^{H} u \, dz \qquad (8) $$

the expression (7) can be written as:

$$ \int_{z_b}^{H} \frac{\partial u}{\partial x} dz = \frac{\partial}{\partial x} [h\bar{u}] - [u]_{H} \left( \frac{\partial h}{\partial x} + \frac{\partial z_b}{\partial x} \right) + [u]_{z_b} \frac{\partial z_b}{\partial x} \qquad (9) $$
• For z-term in (5), the following expression is obtained applying the Leibniz rule:

$$ \int_{z_b}^{H} \frac{\partial w}{\partial z} dz = [w]_{H} - [w]_{z_b} \qquad (10) $$

Using (9) and (10), the expression (5) can be expressed as:

$$ \frac{\partial}{\partial x}(h\bar{u}) - [w]_H \frac{\partial h}{\partial x} - [w]_H \frac{\partial z_b}{\partial x} + [w]_{z_b} \frac{\partial z_b}{\partial x} + [w]_H - [w]_{z_b} = 0 \qquad (11) $$

By using the chain rule to develop the first derived term, $\partial(h\bar{u})/\partial x$, and by using an uniform vertical profile of horizontal velocity as $\bar{u} = [u]_H = [u]_{z_b}$, the expression (11) can be expressed as:

$$ h \frac{\partial \bar{u}}{\partial x} + \bar{u} \frac{\partial h}{\partial x} - \bar{u} \frac{\partial h}{\partial x} - \bar{u} \frac{\partial z_b}{\partial x} + \bar{u} \frac{\partial z_b}{\partial x} + [w]_H - [w]_{z_b} = 0 \qquad (12) $$

and, therefore:

$$ h \frac{\partial \bar{u}}{\partial x} + [w]_H - [w]_{z_b} = 0 \Rightarrow \frac{\partial \bar{u}}{\partial x} + \frac{[w]_H - [w]_{z_b}}{h} = 0 \qquad (13) $$

The main difference between the mass conservation equation in the system of Eqs. (1) and the DIIC, both derived from the divergence-free condition, lies in the assumed velocity profiles and the application of boundary conditions. In the derivation of the DIIC, only the bottom boundary condition, $[w]_{z_b} = \bar{u}(\frac{\partial z_b}{\partial x})$, is applied, while the other boundary conditions are not explicitly enforced. As a result, the vertical velocity remains in the equations. Consequently, the formulation requires not only the definition of a depth-averaged horizontal velocity $\bar{u}$, but also a depth-averaged vertical velocity $\bar{w}$:

$$ \bar{w} = \frac{1}{h} \int_{z_b}^{H} w \, dz \qquad (14) $$

While an uniform profile for horizontal velocity is considered, a linear profile for the vertical velocity along the z-coordinate is taken, with the referenced coordinate $z' = z - z_b$:

$$ w(z') = [w]_{z_b} + \frac{[w]_H - [w]_{z_b}}{h} z' \qquad (15) $$

and inserting Eq. (15) in (14):

$$ \bar{w} = \frac{1}{h} \left[ \int_{0}^{h} w_{z_b} dz' + \int_{0}^{h} \left( \frac{w_H - w_{z_b}}{h} \right) z' dz' \right] \Rightarrow \bar{w} = \frac{1}{h} \left[ [w]_{z_b} h + \left( \frac{[w]_H - [w]_{z_b}}{h} \right) \frac{h^2}{2} \right] \qquad (16) $$

the following expression is obtained:

$$ [w]_H = 2\bar{w} - [w]_{z_b} \qquad (17) $$

Using (17) in (13) and considering the horizontal velocity u uniform along the vertical:

$$ \frac{\partial \bar{u}}{\partial x} + \frac{2\bar{w} - 2[w]_{z_b}}{h} = 0 \qquad (18) $$

Applying the bottom kinematic boundary condition for a static bed, $[w]_{z_b} = \bar{u}(\frac{\partial z_b}{\partial x})$, in (18):

$$ \frac{\partial \bar{u}}{\partial x} + \frac{2\bar{w} - 2\bar{u}(\frac{\partial z_b}{\partial x})}{h} = 0 \Rightarrow \frac{\partial \bar{u}}{\partial x} + \frac{2\bar{w}}{h} - \frac{2}{h}(\bar{u}\frac{\partial z_b}{\partial x}) = 0 \qquad (19) $$

and multiplying by h, (19) can be modified to:

$$ h \frac{\partial \bar{u}}{\partial x} + 2i\bar{w} - 2\bar{u} \frac{\partial z_b}{\partial x} = 0 \qquad (20) $$

Finally, from (20), the expression for the depth-averaged vertical velocity [$LT^{-1}$] is obtained:

$$ \bar{w} = \bar{u} \frac{\partial z_b}{\partial x} - \frac{h}{2} \frac{\partial \bar{u}}{\partial x} \qquad (21) $$

It is important to remark that this vertical velocity is derived exclusively from domain characteristics and conserved variables but is not transported. Therefore, its calculation does not affect the temporal evolution of conservative variables for the 1DH SWE model.

Before incorporating the vertical velocity approximation into the Lagrangian model, its accuracy must be verified. For this reason, and for simplicity, the approximation is developed within a 1DH SWE framework and will be compared with two NHP SWE models in the following sections to assess its validity.

### 3. Numerical scheme for 1DH SWE model

In this section we provide an overview of the numerical scheme for the 1DH SWE model. For more details, the following literature can be used (Morales-Hernández et al., 2013, 2015). For details about the numerical scheme of the 1D NHP models, see Escalante et al. (2019), Echeverribar et al. (2023a). To obtain the updating of the 1DH SWE conservative variables, in the present work a first order explicit upwind numerical scheme is used:

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t^n}{\Delta x} \left[ \sum_{m=1}^{2} (\bar{\lambda}^+ \bar{\gamma} \bar{e})_{i-1/2}^m + \sum_{m=1}^{2} (\bar{\lambda}^- \bar{\gamma} \bar{e})_{i+1/2}^m \right]^n \qquad (22) $$

where $\Delta t = t^{n+1} - t^n$ is the time step [$T$], $\Delta x$ is the cell size [$L$], $\bar{\lambda}$ is the Jacobian matrix eigenvalues, $\bar{e}$ is the Jacobian matrix eigenvectors, $\bar{\gamma}^m_{i+1/2}$ is a compact expression for the wave and source strengths, and $\bar{\lambda}^{\pm m}_{i+1/2} = \frac{1}{2}(\bar{\lambda} \pm |\bar{\lambda}|) \bar{e}^m_{i+1/2}$ accounts for the upwind discretization. The time step $\Delta t^n$ is obtained dynamically and, due to the explicit nature of the numerical scheme, it must be restricted by the Courant–Friedrichs–Lewy (CFL) condition to ensure numerical stability:

$$ \Delta t^n = \text{CFL} \min_{m,k} \left( \frac{\Delta x}{|\bar{\lambda}_k^m|} \right), \quad 0 < \text{CFL} \le 1 \qquad (23) $$

#### 3.1. Details for obtaining the vertical velocity

Considering (21) for the vertical velocity, the spatial derivatives are discretized using a centered method as follows for a general interior cell:

$$ \bar{w}_i^{n+1} = \bar{u}_i^{n+1} \left( \frac{z_{b,i+1} - z_{b,i-1}}{2\Delta x} \right) - \frac{h_i^{n+1}}{2} \left( \frac{\bar{u}_i^{n+1} - \bar{u}_{i-1}^{n+1}}{2\Delta x} \right) \qquad (24) $$

It is important to note that, this expression must be modified at the inlet and outlet boundary nodes. For example, when $i = 0$, (24) is replaced by:

$$ \bar{w}_0^{n+1} = \bar{u}_0^{n+1} \left( \frac{z_{b,1} - z_{b,0}}{\Delta x} \right) - \frac{h_0^{n+1}}{2} \left( \frac{\bar{u}_1^{n+1} - \bar{u}_0^{n+1}}{\Delta x} \right) \qquad (25) $$

As can be seen, when $i$ is equal to 0, the spatial derivatives in x-direction must be solved using a forward method. When $i = N - 1$, the spatial derivatives must be discretized by a backward method:

$$ \bar{w}_{N-1}^{n+1} = \bar{u}_{N-1}^{n+1} \left( \frac{z_{b, N-1} - z_{b, N-2}}{\Delta x} \right) - \frac{h_{N-1}^{n+1}}{2} \left( \frac{\bar{u}_{N-1}^{n+1} - \bar{u}_{N-2}^{n+1}}{\Delta x} \right) \qquad (26) $$

### 4. Analytical cases for the vertical velocity approximation

The following test cases are used to compare the vertical velocity obtained from the proposed 1DH SWE model, using expression (21), with that of the 1D HE-NHP model (27) and the 1D HR-NHP model (28). Since the HE-NHP and HR-NHP models are used to benchmark the proposed approach against more complex formulations, the additional equations required to complete these models are also provided. In the case of the 1D HE-NHP model (Escalante et al., 2018; Bristeau et al., 2015; Yamazaki et al., 2009), the system of equations is:

$$ \left\{ \begin{array}{l} \frac{\partial h}{\partial t} + \frac{\partial(h\bar{u})}{\partial x} = 0 \\ \frac{\partial(h\bar{u})}{\partial t} + \frac{\partial}{\partial x} \left( \frac{h\bar{u})^2}{h} + \frac{1}{2}gh^2 + (hp_{nh}) \right) = -(gh + 2p_{nh}) \frac{\partial z_b}{\partial x} - \frac{\tau_b}{\rho} \\ \frac{\partial(h\bar{u}_{nh})}{\partial t} + \frac{\partial(h\bar{u}\bar{w}_{nh})}{\partial x} = 2p_{nh} \\ h \frac{\partial(h\bar{u})}{\partial x} + 2h\bar{w}\bar{u}_{he-nh} - h\bar{u} \frac{\partial(h+2z_b)}{\partial x} = 0 \end{array} \right. \qquad (27) $$
<figure><img src="image_3.png" /><figcaption>Fig. 1. States at t = 6 s (a), t = 9 s (b), and t = 12 s (c) showing the spatial evolution of the vertical velocity for each model along the longitudinal coordinate. Temporal evolution of the vertical velocity RMSE between the 1DH SWE and 1D HR-NHP models, as well as between the 1DH SWE and 1D HE-NHP models, for the test case 1 (d).</figcaption></figure>

where $p_{nh}$ is the depth-averaged non-hydrostatic pressure correction [$ML^{-1}T^{-2}$], which must be solved implicitly from the last equation using $h$, $\bar{u}$ and $\bar{v}_{nh}$, being $\bar{v}_{nh}$ the vertical velocity obtained by a non-hydrostatic pressure assumption.

For the 1D HR-NHP model (Escalante et al., 2019; Echeverriar et al., 2023a), the system of equations is:

$$ \left\{ \begin{aligned} & \frac{\partial h}{\partial t} + \frac{\partial(h\bar{u})}{\partial x} = 0 \\ & \frac{\partial(h\bar{u})}{\partial t} + \frac{\partial}{\partial x} \left( \frac{(h\bar{u})^2}{h} + \frac{1}{2}gh^2 + (hp_{nh}) \right) = -(gh + 2p_{nh}) \frac{\partial z_h}{\partial x} - \frac{\tau_h}{\rho} \\ & \frac{\partial(h\bar{v}_{nh})}{\partial t} + \frac{\partial(h\bar{u}\bar{v}_{nh})}{\partial x} = 2p_{nh} \\ & \frac{\partial(hp_{nh})}{\partial t} + \frac{\partial}{\partial x} \left( \frac{(hu)(hp_{nh})}{h} \right) + c_p^2 \frac{\partial(hu)}{\partial x} + c_p^2 \left( 2\bar{v}_{nh} - \bar{u} \frac{\partial(h+2z_h)}{\partial x} \right) = 0 \end{aligned} \right. \quad (28) $$

where $c_p = \alpha \sqrt{gH_0}$, being $H_0$ a reference depth [L] and $\alpha$ a parameter that tunes the celerity of the pressure waves (Escalante et al., 2019), where a value of $\alpha = 3$ is used based on the literature (Echeverriar et al., 2023a). Note that, when $c_p \to \infty$, the 1D HE-NHP equations system (27) is recovered. This is because $c_p$ can be interpreted as the pressure waves speed in the depth integrated NHP system (Escalante et al., 2019). For the 1D HE-NHP, the pressure waves propagation is considered instantaneous.

### 4.1. Test case 1: Dam-break

The first test case consists of a flat and frictionless horizontal channel with a constant width of $B = 20$ m and a length of $L = 300$ m, with initial discontinuous water depth and zero velocity. The water depth is 4 m on the left region and 1 m on the right region of an ideal gate that opens suddenly at $t = 0$. The cell size for this case is $\Delta x = 0.25$ m. The goal is to simulate a dam-break scenario, highlighting the vertical velocity variations in such situations. The spatial evolution of the vertical velocity at different states is shown in Fig. 1. As observed in this figure, the differences between the approximation of the 1DH SWE model and the 1D NHP models are generally insignificant, except at the highest vertical velocity peak. Moreover, the propagation behavior of the models is very similar in both the left and right directions.

To quantify the error between both models, the Root Mean Squared Error (RMSE) is computed as (Barnston, 1992):

$$ \text{RMSE}(t^n) = \sqrt{\frac{1}{N_L} \sum_{l=1}^{N_L} (\bar{v}_l(t^n) - \bar{v}_{nh,l}(t^n))^2} \quad (29) $$

where $N_L$ is the number of measurement points along the longitudinal coordinate, $\bar{v}_l(t^n)$ is the vertical velocity of the 1DH SWE model at the
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** (a) Initial condition for the reef test case. (b-e) Spatial profiles of vertical velocity at $t = 1.5$ s, 3.0 s, 4.5 s, and 6.0 s, respectively, comparing the 1DH SWE model with vertical velocity and the 1D NHP models along the channel.</figcaption></figure>

spatial position $l$ and time $t^n$, and $\bar{w}_{nh,l}(t^n)$ is the vertical velocity of the 1D NHP models at spatial position $l$ and time $t^n$. The temporal evolution of the vertical velocity RMSE between the 1DH SWE model and the 1D NHP models is shown in Fig. 1(d). As seen in this figure, the RMSE decreases over time and stabilizes when the vertical velocity peaks no longer change.

## 4.2. Test case 2: Reef

The second case for the vertical velocity accuracy assessment is conducted using a well-known test case described in Escalante et al. (2018) (see Fig. 2(a)). The setup consists of a channel with a length of $L = 45$ m and a constant width of $B = 9$ m. The channel roughness is uniform, with a Manning's coefficient $n = 0.01 \text{ m}^{-1/3}$. The upstream boundary condition is a constant inlet discharge $Q_{in} = 1 \text{ m}^3 \text{ s}^{-1}$, imposed at the left boundary wall, while the remaining boundary walls impose a reflective boundary condition. The initial condition is

a constant water surface level of 1 m. The cell size for this case is $\Delta x = 0.045$ m. The spatial evolution of the vertical velocity at different times is shown in Fig. 2. As seen in this figure, the differences between the 1DH SWE model and the 1D NHP models are significant at the beginning of the simulation. However, as the flow reaches the dry area, these differences diminish. Moreover, the propagation behavior of the models is very similar both in the left and right directions. The temporal evolution of the vertical velocity RMSE is shown in Fig. 3a, where the RMSE value progressively decreases as the flow propagates. This trend is further illustrated in Fig. 3b, which presents the maximum vertical velocity for each model, showing that the values converge and become really similar by the end of the simulation.

## 4.3. Test case 3: Plane beach

The third case for assessing vertical velocity accuracy is conducted using another well-known test case (see Fig. 4(a)). The details of
<figure><img src="image_5.png" /><figcaption>Fig. 3. Temporal evolution of the vertical velocity RMSE between the 1DH SWE and 1D HR-NHP models, as well as between the 1DH SWE and 1D HE-NHP models, for the reef test case (a). Temporal evolution of the maximum vertical velocity value provided by the different models for the test case 2 (b).</figcaption></figure>

the analytical setup are described in Escalante et al. (2018), which considers a channel with a length of L = 49.85 m and a constant width of B = 10 m. The channel roughness is uniform, with a Manning's coefficient n = 0.01 s m⁻¹/³. The upstream boundary condition consists of a constant inlet discharge Qᵢₘ = 3 m³ s⁻¹, imposed at the left boundary wall, while at the outlet boundary a reflective boundary condition is imposed. The initial condition is a constant water surface level of 1 m. The cell size for this case is Δx = 0.05 m. Considering that the results of the NHP models are very similar (see Figs. 1–3), in this case, only the 1DH SWE model is compared with the 1D HR-NHP model. Some selected temporal states and the temporal evolution of the RMSE for vertical velocity are shown in Fig. 4. As observed in the two previous test cases, a significant difference between the models appears in the peak vertical velocity during the initial moments (see Fig. 4(b)). However, as the peak loses amplitude and the vertical velocity becomes less relevant to the flow evolution, this difference decreases until negligible errors are obtained, as shown in Figs. 4(d) and 4(e).

## 5. Extension of the governing equations to two-dimensions

After verifying the accuracy of the vertical velocity approximation in one dimension, the following section extends the governing equations to two dimensions to improve the Lagrangian model for particle transport. Specifically, it presents the 2DH SWE governing the flow, along with the equations defining the Lagrangian particle transport.

### 5.1. 2DH SWE model

The two-dimensional hydrodynamic model describes the depth-averaged surface flow by means of the hyperbolic 2DH SWE system, is based on mass and momentum conservation under hydrostatic pressure hypothesis:

$$ \frac{\partial U}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y} = S_b + S_f \quad (30) $$

where:

$$ \mathbf{U} = \begin{pmatrix} h \\ h\bar{u} \\ h\bar{v} \end{pmatrix}, \quad \mathbf{F} = \begin{pmatrix} h\bar{u} \\ \bar{u}^2 + \frac{1}{2}gh^2 \\ h\bar{u}\bar{v} \end{pmatrix}, \quad \mathbf{G} = \begin{pmatrix} h\bar{v} \\ h\bar{u}\bar{v} \\ \bar{v}^2 + \frac{1}{2}gh^2 \end{pmatrix}, \quad (31) $$

$$ S_b = \begin{pmatrix} 0 \\ -gh \frac{\partial z_b}{\partial x} \\ -gh \frac{\partial z_b}{\partial y} \end{pmatrix}, \quad S_f = \begin{pmatrix} 0 \\ -\frac{gh n^2 \bar{u} \sqrt{\bar{u}^2 + \bar{v}^2}}{h^{4/3}} \\ -\frac{gh n^2 \bar{v} \sqrt{\bar{u}^2 + \bar{v}^2}}{h^{4/3}} \end{pmatrix} \quad (31) $$

being h the water depth [L], u = (ā, ā, 0) the velocity flow vector [LT⁻¹], g the gravitational acceleration [LT⁻²], z_b the bottom level [L], and n the Manning's roughness coefficient [TL⁻¹/³] (Arcement and Schneider, 1984).

### 5.2. Particle tracking equations in 2DH problems

For particles with a negligible mass/volume, transport is driven by advection and diffusion, being the evolution of the particle position x_p = (x_p, y_p, z_p) governed by the horizontal velocity components:

$$ \left\{ \begin{aligned} \frac{dx_p}{dt} &= \bar{u}(x_p) + u_{\text{disp}}(x_p) \\ \frac{dy_p}{dt} &= \bar{v}(x_p) + v_{\text{disp}}(x_p) \\ \frac{dz_p}{dt} &= 0 \end{aligned} \right. \quad (32) $$

where v_{disp} = (u_{disp}, v_{disp}, 0) is the horizontal velocity induced by the dispersion. This velocity can be modeled using various dispersion models. In this study, a random-walk model is employed to simulate dispersion (Jalón-Rojas et al., 2019). This model proposes to formulate v_{disp} as a function of the diffusivity coefficients K_{hx} and K_{hy} [L²T⁻¹] in the x- and y-coordinates, respectively. In contrast to imposing constant values based on empirical results (Jalón-Rojas et al., 2019; Peeters and Hofmann, 2015) that depend on the size and shape of the object, the diffusivity coefficients K_{hx} and K_{hy} are here derived from a velocity-dependent expression, as turbulence is directly proportional to velocity. Therefore, the standard anisotropic diffusion model (Vallés et al., 2023) is utilized to determine the diffusivity variables, with the following expressions:

$$ K_{hx} = \epsilon_L |u^*| h(x_p), \quad K_{hy} = \epsilon_T |u^*| h(x_p) \quad (33) $$

where ε_L and ε_T are the longitudinal and transversal coefficients (Rutherford, 1994), respectively; and |u^*| is the friction velocity [LT⁻¹]:

$$ |u^*| = \sqrt{\frac{\tau_o}{\rho}} \quad (34) $$

being τ_o the bottom shear stress [ML⁻¹T⁻²] and ρ the water density [ML⁻³]. As can be seen in (31), the bottom shear stress is modeled by the Manning's roughness law. Therefore, the friction velocity is defined by:

$$ |u^*| = n \sqrt{g \frac{\bar{u}^2 + \bar{v}^2}{(h(x_p))^{1/3}}} \quad (35) $$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Initial condition for the plane beach test case (a). States at t = 3.0 s (b), t = 4.5 s (c), and t = 6.0 s (d), showing the spatial evolution of the vertical velocity for each model along the longitudinal coordinate. Temporal evolution of the vertical velocity RMSE between the 1DH SWE and 1D HR-NHP models for the test case 3 (e).</figcaption></figure>

### 5.3. Adding the vertical velocity in the 2DH SWE model

For the 2DH SWE model, the system of Eqs. (30) is formed by the depth integrated mass conservation and two momentum equations in the x- and y-directions. The integration of the divergence-free condition between the bottom and the water surface level is defined as:

$$ \int_{z_b}^{H} \frac{\partial u}{\partial x} dz + \int_{z_b}^{H} \frac{\partial v}{\partial y} dz + \int_{z_b}^{H} \frac{\partial w}{\partial z} dz = 0 \quad (36) $$

where (u, v, w) are the velocity components [LT⁻¹] in the x-, y- and z-coordinates, respectively. Using the Leibniz rule, the terms of the expression (36) can be expressed as:

* For the x-term in (36), the following expression is obtained applying the Leibniz rule as is detailed in Section 2.2:

$$ \int_{z_b}^{H} \frac{\partial u}{\partial x} dz = \frac{\partial}{\partial x}(h\bar{u}) - [u]_H \left( \frac{\partial h}{\partial x} + \frac{\partial z_b}{\partial x} \right) + [u]_{z_b} \frac{\partial z_b}{\partial x} \quad (37) $$

* For the y-term in (36), the following expression is obtained applying the Leibniz rule and operating analogous to the x-term in Section 2.2:

$$ \int_{z_b}^{H} \frac{\partial v}{\partial y} dz = \frac{\partial}{\partial y}(h\bar{v}) - [v]_H \left( \frac{\partial h}{\partial y} + \frac{\partial z_b}{\partial y} \right) + [v]_{z_b} \frac{\partial z_b}{\partial y} \quad (38) $$
• For the z-term, the following expression is obtained applying the Leibniz rule:

$$\int_{z_b}^{H} \frac{\partial w}{\partial z} dz = [w]_H - [w]_{z_b} \qquad (39)$$

Using (37), (38) and (39), the expression (36) can be expressed as:

$$\frac{\partial}{\partial x}(h\bar{u}) - [u]_H \left(\frac{\partial h}{\partial x} + \frac{\partial z_b}{\partial x}\right) + [u]_{z_b} \frac{\partial z_b}{\partial x} + \frac{\partial}{\partial y}(h\bar{v}) - [v]_H \left(\frac{\partial h}{\partial y} + \frac{\partial z_b}{\partial y}\right) + [v]_{z_b} \frac{\partial z_b}{\partial y} + [w]_H - [w]_{z_b} = 0 \qquad (40)$$

By using the chain rule to develop the derived terms $\partial(h\bar{u})/\partial x$ and $\partial(h\bar{v})/\partial y$, and by using an uniform vertical profile of the horizontal and transversal velocities as $\bar{u} = [u]_H = [u]_{z_b}$ and $\bar{v} = [v]_H = [v]_{z_b}$, the expression (40) can be expressed as:

$$h \frac{\partial \bar{u}}{\partial x} + h \frac{\partial \bar{v}}{\partial y} + [w]_H - [w]_{z_b} = 0 \Rightarrow \frac{\partial \bar{u}}{\partial x} + \frac{\partial \bar{v}}{\partial y} + \frac{[w]_H - [w]_{z_b}}{h} = 0 \qquad (41)$$

Defining a linear evolution for the vertical velocity $w$ along the z-coordinate as is detailed in Section 2.2:

$$w(z') = [w]_{z_b} + \frac{[w]_{h_s} - [w]_{z_b}}{h} z' \qquad (42)$$

where $z' \in [0, h]$ is the vertical position from the bottom elevation, the following expression can be obtained:

$$[w]_H = 2\bar{w} - [w]_{z_b} \qquad (43)$$

and using (43), the expression (41) becomes:

$$\frac{\partial \bar{u}}{\partial x} + \frac{\partial \bar{v}}{\partial y} + \frac{2\bar{w} - 2[w]_{z_b}}{h} = 0 \qquad (44)$$

Applying the bottom boundary condition, $[w]_{z_b} = \bar{u}(\frac{\partial z_b}{\partial x}) + \bar{v}(\frac{\partial z_b}{\partial y})$, in (44):

$$\frac{\partial \bar{u}}{\partial x} + \frac{\partial \bar{v}}{\partial y} + \frac{2\bar{w}}{h} - \frac{2}{h} \left( \bar{u} \frac{\partial z_b}{\partial x} + \bar{v} \frac{\partial z_b}{\partial y} \right) = 0 \qquad (45)$$

Finally, the expression for the vertical velocity component $\bar{w}$ is obtained:

$$\bar{w} = \bar{u} \frac{\partial z_b}{\partial x} + \bar{v} \frac{\partial z_b}{\partial y} - \frac{h}{2} \frac{\partial \bar{u}}{\partial x} - \frac{h}{2} \frac{\partial \bar{v}}{\partial y} \qquad (46)$$

This expression can be obtained by another procedure, as shown in Ferrari and Saleri (2004).

## 5.4. Improving vertical velocity expression

From the linearity assumption of the vertical component, the vertical velocity at the free surface is defined by Eq. (43). Assuming a null vertical velocity at the bottom ($[w]_{z_b} = 0$), the expression for the vertical velocity $w(z)$ is:

$$w(z) = \frac{2\bar{w}}{h}(z - z_b) \quad \text{if } z_b \le z \le h_s \qquad (47)$$

This is useful when a particle is transported by the flow, having different positions in the vertical, and therefore, being transported by a non-uniform velocity field in this coordinate.

## 6. Numerical schemes for 2D models

In this section we provide an overview of the numerical schemes implemented in SERGHEI for the SWE module (for more details, see Caviedes-Voullière et al. (2023)) and for the LPT module (for more details, see Vallés et al. (2023, 2025)).

### 6.1. Numerical scheme for 2DH SWE model

SERGHEI-SWE uses a first-order accurate upwind finite-volume scheme with a forward Euler time integration to solve the system (30)

on uniform Cartesian grids with grid spacing $\Delta x$ [$L$]. Well-balancing and water depth positivity are ensured by solving numerical fluxes at each cell edge $k$ with augmented Riemann solvers based on the Roe linearization. In fluctuation form, the rule for updating the conserved variables in cell $i$ from time step $n$ to time step $n+1$ reads as follows:

$$\mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t^n}{\Delta x} \sum_{k=1}^{4} \sum_{m=1}^{3} \frac{\tilde{\lambda}_k^{-}}{\tilde{\lambda}_k} [\tilde{\lambda}_k \tilde{\alpha}_k - \tilde{\beta}_k] \bar{e}_{m,k}^n \qquad (48)$$

where $\tilde{\lambda}_k$ and $\bar{e}_k$ are the eigenvalues and eigenvectors of the linearized system of equations, $\tilde{\alpha}_k$ and $\tilde{\beta}_k$ are the fluxes and bed slope and friction source term linearizations, respectively, and the minus sign for $\tilde{\lambda}_k$ accounts for the upwind discretization.

The time step $\Delta t^n$ is obtained dynamically and restricted by the Courant–Friedrichs–Lewy (CFL) condition to ensure numerical stability:

$$\Delta t^n = \text{CFL} \min_i \left( \frac{\Delta x}{\left| \frac{\tilde{h}_i}{h} \right|^n + \sqrt{g \tilde{h}_i^n}, \left| \frac{\tilde{h}_i}{h} \right|^n + \sqrt{g \tilde{h}_i^n}} \right), \quad 0 < \text{CFL} \le 0.5 \qquad (49)$$

### 6.1.1. Details for obtaining the vertical velocity

The discretization of (46) is performed again using a centered approach. Two indexes $(i, j)$ will be used in this subsection to show how the spatial derivatives are calculated according to the x- or y-direction. These indexes have the ranges of values $i = 0, \dots, N-1$, being $N$ the number of the mesh columns, and $j = 0, \dots, M-1$, being $M$ the number of the mesh rows. Using this notation, the expression (46) is discretized as follows:

$$\bar{w}_{(i,j)}^{n+1} = \bar{u}_{(i,j)}^{n+1} \left( \frac{z_b(i+1,j) - z_b(i-1,j)}{2\Delta x} \right) + \bar{v}_{(i,j)}^{n+1} \left( \frac{z_b(i,j+1) - z_b(i,j-1)}{2\Delta x} \right) - \frac{h_{(i,j)}^{n+1}}{2} \left( \frac{\bar{u}_{(i+1,j)}^{n+1} - \bar{u}_{(i-1,j)}^{n+1}}{2\Delta x} \right) - \frac{h_{(i,j)}^{n+1}}{2} \left( \frac{\bar{v}_{(i,j+1)}^{n+1} - \bar{v}_{(i,j-1)}^{n+1}}{2\Delta x} \right) \qquad (50)$$

However, this expression must be modified when the indexes have an extreme value $(0, N-1, \text{ or } M-1)$. For example, when $i = 0$, (50) is modified by:

$$\bar{w}_{(0,j)}^{n+1} = \bar{u}_{(0,j)}^{n+1} \left( \frac{z_b(i,j) - z_b(0,j)}{\Delta x} \right) + \bar{v}_{(0,j)}^{n+1} \left( \frac{z_b(0,j+1) - z_b(0,j-1)}{2\Delta x} \right) - \frac{h_{(0,j)}^{n+1}}{2} \left( \frac{\bar{u}_{(1,j)}^{n+1} - \bar{u}_{(0,j)}^{n+1}}{\Delta x} \right) - \frac{h_{(0,j)}^{n+1}}{2} \left( \frac{\bar{v}_{(0,j+1)}^{n+1} - \bar{v}_{(0,j-1)}^{n+1}}{2\Delta x} \right) \qquad (51)$$

When $i$ is equal to 0, the spatial derivatives in x-direction must be solved using an forward method. The same modification applies for the y-direction when $j = 0$.

When $i = N-1$, the spatial derivatives in x-direction of the expression (50) must be solved by a backward method, being the same strategy when $j = M-1$:

$$\bar{w}_{(N-1,j)}^{n+1} = \bar{u}_{(N-1,j)}^{n+1} \left( \frac{z_b(N-1,j) - z_b(N-2,j)}{\Delta x} \right) + \bar{v}_{(N-1,j)}^{n+1} \left( \frac{z_b(N-1,j+1) - z_b(N-1,j-1)}{2\Delta x} \right) - \frac{h_{(N-1,j)}^{n+1}}{2} \left( \frac{\bar{u}_{(N-1,j)}^{n+1} - \bar{u}_{(N-2,j)}^{n+1}}{\Delta x} \right) - \frac{h_{(N-1,j)}^{n+1}}{2} \left( \frac{\bar{v}_{(N-1,j+1)}^{n+1} - \bar{v}_{(N-1,j-1)}^{n+1}}{2\Delta x} \right) \qquad (52)$$

### 6.2. Numerical scheme for LPT module

The temporal evolution of particle positions is obtained by considering the particles as mathematical points denoted by $\mathbf{p}$, which are advected in space $\mathbf{x}$ following a velocity field $\mathbf{v}$. Although the flow still has no turbulence model, the dispersive effects of turbulence are modeled on the transport of particles is included understanding that the model might lead to unrealistic particle transport if this phenomenon
<figure><img src="image_7.png" /><figcaption>Fig. 5. Particle transport algorithm based on particle position in an example case, where red trajectories indicate incorrect transport, while green trajectories represent a possible transport path.</figcaption></figure>

is not considered (Merritt and Wohl, 2002). These effects are added using a random-walk model, which can generate sufficient dispersion to introduce turbulent motion (Rutherford, 1994):

$$ 
\left\{ 
\begin{aligned}
x_p^{n+1} &= x_p^n + \Delta t_p^n \bar{u}_i^{n+1} + R_{p,x}^n (2\sigma^{-1} K_{h,x,i}^n \Delta t_p^n)^{1/2} \\
y_p^{n+1} &= y_p^n + \Delta t_p^n \bar{e}_i^{n+1} + R_{p,y}^n (2\sigma^{-1} K_{h,y,i}^n \Delta t_p^n)^{1/2} \\
z_p^{n+1} &= z_p^n + \Delta t_p^n \omega (z_p^n)^{n+1} + R_{p,z}^n (2\sigma^{-1} K_{v,i}^n \Delta t_p^n)^{1/2}
\end{aligned}
\right. 
\quad (53) $$

where $v_i^n$ is the Eulerian velocity field in cell $i$ at time step $n$, in which the particle $p$ is located, i.e., contains the point $x_p^n = (x_p^n, y_p^n, z_p^n)$; $h_i^n$ is the water depth in the cell $i$ at time step $n$, $z_{b,i}$ is the bottom elevation of the cell $i$, and $\Delta t_p^n$ is the particle time step. $R_{p,x}^n$, $R_{p,y}^n$ and $R_{p,z}^n$ are random numbers, which follow a uniform distribution with mean 0 and standard deviation $\sigma = 1$. $K_{h,x,i}^n$ and $K_{h,y,i}^n$ are the horizontal diffusivity in $\hat{x}$ and $\hat{y}$ components, respectively, and $K_{v,i}^n$ is the vertical diffusivity. To ensure a proper assessment of the vertical velocity approaches, turbulence effects are neglected by setting the Manning's roughness coefficient to zero in the test cases.

When the vertical position of the particle is likely to be variable, incorporating vertical velocity into the particle trajectories is essential for ensuring a physically consistent transport. As observed in Fig. 5, the predicted particle position using the updating scheme in (53) is altered in the following cases:

1. When the updated horizontal coordinates correspond to a solid body position (condition 2 in Fig. 5), the horizontal velocity updating is not applied, and vertical velocity dominates the particle transport:

$$ x_p^{n+1} = x_p^n, \quad y_p^{n+1} = y_p^n \quad (54) $$

2. When the updated vertical coordinate is above the water column (condition 3 in Fig. 5), the vertical velocity is not applied

because it would result in an unrealistic particle position:

$$ z_p^{n+1} = h_i + z_{b,i} \quad (55) $$

3. When the updated vertical coordinate is below the bed level (condition 4 in Fig. 5), the vertical velocity is not applied because it would result in an unrealistic particle position and the particle is deposited:

$$ z_p^{n+1} = z_{b,i} \quad (56) $$

Therefore, the vertical position of the particle is limited by:

$$ 
\left\{ 
\begin{aligned}
z_p^{n+1} &= h_i^{n+1} + z_{b,i}, & \text{if } z_p^n + \Delta t_p^n \omega (z_p^n)^{n+1} > h_i^{n+1} + z_{b,i} \\
z_p^{n+1} &= z_{b,i}, & \text{if } z_p^n + \Delta t_p^n \omega (z_p^n)^{n+1} < z_{b,i}
\end{aligned}
\right. 
\quad (57) $$

## 7. Test cases for LPT module

The validation of the vertical velocity is conducted through the simulation of several test cases, assessing both accuracy and computational efficiency. The following cases present different non-trivial geometries for particle transport in a vertically averaged 2D model. Thus, they illustrate how particles trajectories are improved with the vertical velocity model. The simulations are run using a GPU NVIDIA GeForce RTX 3060 Ti. To enable a fair comparison between the two vertical velocity approximations, the test cases are configured without bed roughness, ensuring that the dispersion terms of the LPT model, which are inherently stochastic, remain inactive.

### 7.1. Test case 4: bump obstacle

The following case is a flat, frictionless, rectangular channel ($L = 25$ m, $B = 0.1$ m) with a symmetric and smooth bump centered at $x = 10$ m across the whole channel, described in more detail in Delestre et al. (2013). The initial condition is a uniform still water surface level of
<figure><img src="image_8.png" /><figcaption>Fig. 6. Trajectories and selected temporal states of particle transport using the average vertical velocity (a) and the vertically dependent velocity $w(z)$ (b) for the test case 4.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. Temporal evolution of the mean vertical position difference between trajectories computed with the average vertical velocity and those using $w(z)$ for three different initial vertical position ranges (a), and particle trajectories and selected temporal states of particle transport without the vertical velocity approximation (b) for the test case 4.</figcaption></figure>

1 m. The upstream boundary condition is a constant inlet discharge $Q_{in} = 0.110 \text{ m}^3\text{s}^{-1}$, while the downstream boundary condition maintains a constant water surface level of 1 m. The remaining boundaries are solid walls. The initial particle distribution consists of a vertical column of particles (see Fig. 6), spaced equidistantly at 0.05 m intervals. As there is no roughness in the channel, the dispersion terms in particle transport are null. Particle transport is simulated using both the average vertical velocity expression (46) (see Fig. 6a) and the vertical velocity profile $w(z)$ given by (47) (see Fig. 6). These figures show the particles trajectories and four specific temporal states for both vertical velocity approximations. As observed, the particles trajectories at the highest and lowest vertical positions differ between the two velocity expressions. In Fig. 6a, particles at the lowest positions suffer a delay compared to the rest because their horizontal velocity is not added when their updated position falls below the bump height. However, as shown in Fig. 6b, these particles are deposited when the vertical velocity $w(z)$ is applied, as their vertical position is sufficiently close to the bottom elevation, leading to a null vertical velocity. Moreover, vertical transport is more pronounced for particles using $w(z)$ when their vertical position exceeds half the water depth, resulting in a greater difference between the free surface and their final vertical positions. This behavior is also illustrated in Fig. 7a, which shows the temporal evolution of the mean vertical position difference between both trajectories, divided into three different initial vertical position ranges. As shown in the figure, the difference remains negligible in the central region ($h/3 \le z_0 < 2h/3$) throughout the entire simulation. However, in the upper and lower regions, noticeable differences arise

due to the disparity between the average vertical velocity and the vertically dependent velocity profile, $w(z)$. When the particles pass over the bump obstacle, the differences follow a similar pattern: in the upper region, particles using $w(z)$ exhibit a greater upward displacement than those using the average vertical velocity, while in the lower region, they exhibit a smaller downward displacement. The vertical velocity approximation in both versions improves particle trajectories, as shown in Fig. 7b, which displays trajectories without any vertical velocity correction. As observed, the trajectories are influenced solely by horizontal advection; consequently, particles with vertical positions lower than the bump height are not able to cross it. Moreover, particles with higher vertical positions are not affected by the presence of the bump and do not experience any displacement.

## 7.2. Test case 5: Sharp bottom obstacle with a cavity

The second test case for the LPT model consists of a flat, frictionless, rectangular channel ($L = 1 \text{ m}, B = 0.15 \text{ m}$) with a sharp obstacle at the bottom across the section that contains a vertical cavity at its midpoint, similar to the experiment presented in You and Tinoco (2023). As sketched in Fig. 8a, the obstacle is 0.105 m in height and the cavity is 0.04 m in the longitudinal direction. The initial condition is uniform still water surface level of 1.3 m. The upstream boundary condition is a constant inlet discharge $Q_{in} = 0.442 \text{ m}^3\text{s}^{-1}$, while the downstream boundary condition is a constant water surface level of 1.3 m. All remaining walls impose a reflective boundary condition. As there is no roughness along the channel, the dispersion terms in
<figure><img src="image_10.png" /><figcaption>Fig. 8. Initial flow condition (a) and the average vertical velocity and bottom elevation along the longitudinal coordinate at the center of the channel (b) for the test case 5.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Particle trajectories and selected temporal states using the average vertical velocity (a) and the vertically dependent velocity $w(z)$ (b) for the test case 5.</figcaption></figure>

particle transport are null. The initial particle distribution consists of a vertical column of particles (see Fig. 9), spaced evenly at 0.01 m intervals. The flow reaches a steady state, with the spatial evolution of the average vertical velocity and the water depth shown in Fig. 8b. As illustrated, the vertical velocity remains approximately zero along the channel, except near the boundaries of the vertical cavity walls. The water surface level and particles trajectories computed using the average vertical velocity are shown in Fig. 9a, while those using the vertically dependent velocity $w(z)$ are shown in Fig. 9b. The differences between the trajectories in Figs. 9a and 9b resemble those observed in the first test case. However, in this scenario, the details of the algorithm are demonstrated more clearly. When the average velocity is used, particles undergo the same ascending or descending motion regardless of their vertical position, which leads to the deposition of some particles at the bottom of the first vertical cavity wall (similar to case 4 in Fig. 5). In contrast, when the vertically dependent velocity is used, particles take longer to pass over the vertical cavity wall, experiencing a delay compared to particles with an initial vertical position above the height of the vertical cavity walls. These particles are not deposited, as their vertical positions are nearly equal to the vertical cavity wall bottom, and therefore, the vertical velocity $w(z)$ is negligible. Moreover, as shown in test case 4, the vertical velocity approximation in both versions improves particle trajectories, as illustrated in Fig. 10, which again displays trajectories without any vertical velocity correction. In this case, the trajectories result solely from horizontal advection. Consequently, as in the previous test case, particles with vertical positions lower than the cavity wall height are

not able to cross it. Furthermore, no particles are transported into the cavity, which differs from the behavior reported in You and Tinoco (2023).

### 7.3. Test case 6: Conical island

The next test case involves unsteady flow around a conical island, as presented in Briggs et al. (1995). The geometry consists of a truncated cone with a base diameter of 7.2 m, a top diameter of 2.2 m, and a height of 0.625 m, located at the center of a smooth, flat domain measuring 25 m × 30 m (see Fig. 11). A uniform Manning roughness coefficient of $n = 0.013 \text{ s m}^{-1/3}$ is imposed throughout the domain. The initial condition is still uniform water level of $h_0 = 0.32 \text{ m}$, and a wave is imposed at the inlet boundary side following the expression for the water depth:

$$h(t) = h_0 + A \operatorname{sech}^2 \left( \frac{B(t-T)}{C} \right), \quad B = \sqrt{gh_0} \left( 1 + \frac{A}{2h_0} \right), \\ C = h_0 \sqrt{\frac{4h_0 B}{3A \sqrt{gh_0}}} \quad (58)$$

where $A = 0.042 \text{ m}$ is the wave amplitude, and $T = 2.84 \text{ s}$ is the time at which peak of the wave enters the domain. The remaining boundaries are treated as solid walls. The initial condition for the 250 particles is a random distribution all over the wet domain. The temporal evolution of particle trajectories for both vertical velocity approximations are shown in Figs. 12 and 13. As observed, the trajectories of some particles differ
<figure><img src="image_12.png" /><figcaption>Fig. 10. Particle trajectories and selected temporal states of particle transport without the vertical velocity approximation for the test case 5.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Initial flow condition for the test case 6.</figcaption></figure>

significantly between the two approaches. This is because, when the average vertical velocity is used, some particles are deposited and come to rest. These differences are further illustrated in Fig. 14, which shows the temporal evolution of the mean vertical position difference between both approximations, categorized into three ranges based on the initial vertical position of each particle. The behavior of these differences is consistent with the trends observed in the previous test cases.

## 8. Conclusions

This work presents an approximation for the vertical flow velocity aimed at enhancing a Lagrangian particle transport in 2DH SWE hydraulic models, without requiring the high computational cost associated with fully 3D models.

First, three one-dimensional test cases are analyzed to assess the accuracy of the vertical velocity approximation, comparing results with those from two NHP models. In all cases, the differences between the three models are small and become negligible once the flow reaches steady state. However, in transient situations, such as the initial moments of the test case 1 and 3, some discrepancies are observed, particularly in the peak vertical velocities. These differences arise because NHP models account for the dynamic coupling between vertical velocity and pressure, which influences the evolution of the conserved variables. In contrast, the proposed approximation does not capture this instantaneous interaction. Nevertheless, the results confirm that the vertical velocity approximation can accurately represent vertical flow behavior without solving the more complex equations required by NHP models or incurring the additional computational cost of explicitly transporting the vertical velocity.

Next, the vertical velocity approximation is applied to particle transport simulations. The average vertical velocity is improved assuming

<figure><img src="image_14.png" /><figcaption>Fig. 12. X-Y plane trajectories and initial particle state using the average vertical velocity: full domain view (a) and zoomed-in view of the central area (b) for the test case 6.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 13. X-Y plane trajectories and initial particle state using the vertically dependent velocity $w(z)$: full domain view (a) and zoomed-in view of the central area (b) for the test case 6.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Temporal evolution of the mean vertical position difference between trajectories computed with the average vertical velocity and those using $w(z)$, for three different initial vertical position ranges in the test case 6.</figcaption></figure>

a vertical linear profile, which proves beneficial for modeling particle movement. To evaluate both vertical velocity approaches, three test cases are simulated to compare particle transport behavior under both the average and vertically dependent velocity approximations. While both methods yielded similar trajectories in the central region of the water column, notable differences are observed in the upper and lower regions of the water column. Specifically, the average velocity approximation causes an excessive particle deposition, whereas the vertically dependent velocity $w(z)$ produces more realistic particle trajectories. In contrast, simulations that do not use any vertical velocity approximation result in particle paths governed solely by horizontal advection. In such cases, particles with vertical positions below the obstacle height are unable to cross it, while those above remain unaffected by the presence of obstacles (see Figs. 7b and 10). Both the average and vertically dependent approximations allow particles to traverse obstacles

(see Figs. 6 and 9); however, the use of $w(z)$ results in more coherent behavior, avoiding unrealistic accumulation.

In summary, the proposed vertical flow approximations are sufficiently accurate to capture the effects of vertical motion without requiring its transport or the use of more complex governing equations, such as those found in NHP models. This allows for a more efficient yet physically meaningful representation of vertical dynamics. Incorporating either the depth-averaged vertical velocity or the vertically varying profile $w(z)$ significantly improves the realism of particle trajectories, with the latter offering superior consistency across test cases. Moreover, the method establishes a practical foundation for applications such as Lagrangian sediment transport modeling, where vertical flow components can play a critical role, all while maintaining compatibility with standard shallow water frameworks.
## CRediT authorship contribution statement

**P. Vallés:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **J. Segovia-Burillo:** Writing – review & editing, Writing – original draft, Visualization, Software. **M. Morales-Hernández:** Writing – review & editing, Writing – original draft, Supervision, Software, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **V. Roeber:** Writing – review & editing, Writing – original draft, Supervision, Methodology, Investigation, Formal analysis, Conceptualization. **P. García-Navarro:** Writing – review & editing, Writing – original draft, Supervision, Resources, Project administration, 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: Pablo Valles reports financial support was provided by University of Zaragoza. Pablo Valles reports financial support was provided by Campus Iberus Consortium. Pablo Valles, Mario Morales-Hernandez, Jose Segovia-Burillo, Pilar Garcia-Navarro reports financial support was provided by Spain Ministry of Science and Innovation. Pablo Valles, Mario Morales-Hernandez, Jose Segovia-Burillo, Pilar Garcia-Navarro reports financial support was provided by Government of Aragón. Mario Morales-Hernandez, Pilar Garcia-Navarro reports financial support was provided by University of Zaragoza. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgments

Pablo Vallés is funded by the UPPA-UNIZAR Research Grant PRD/2022-03 and he was funded by ERASMUS+ KA103 "IBERUS+" 2021-1-ES01-KA130-HED-000004265. This work was supported by project PID2022-137334NB-I00 funded by MCIN/AEI/10.13039/501100011033 and by ERDF/EU. This work has been partially funded by the Government of Aragón, through the research grant T32_23R Fluid Dynamics Technologies. This work was partially funded by the project JIUZ2023-IA-04 by UNIZAR.

## Data availability

Data will be made available on request.

## References

Arcement, G., Schneider, V., 1984. Guide for selecting manning's roughness coefficients for natural channels and flood plains. (2339) U.S. Geol. Surv. Water-Supply Paper.

Armanini, A., Fraccarollo, L., Rosatti, G., 2009. Two-dimensional simulation of debris flows in erodible channels, computers & geosciences. Model. Simul. Danger. Phenom. Hazard Mapp. 35 (5), 993–1006. http://dx.doi.org/10.1016/j.cageo.2007.11.008.

Baharvand, S., Ahmari, H., Taghvaei, P., 2023. Developing a lagrangian sediment transport model for open channel flows. Int. J. Sediment Res. 38 (2), 153–165. http://dx.doi.org/10.1016/j.ijsrc.2022.09.003.

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.

Briggs, M.J., Synolakis, C.E., Harkins, G.S., Green, D.R., 1995. Laboratory experiments of tsunami runup on a circular island. Pure Appl. Geophys. 144 (3), 569–593. http://dx.doi.org/10.1007/BF00874384.

Bristeau, M.-O., Mangeney, A., Sainte-Marie, J., Seguin, N., 2015. An energy-consistent depth-averaged euler system: Derivation and properties. Discret. Contin. Dyn. Syst. - B 20 (4), 961–988.

Burton, R.R., Woodhouse, M.J., Gadian, A.M., Mobbs, S.D., (2020). The use of a numerical weather prediction model to simulate near-field volcanic plumes. Atmosphere. 11. 6. http://dx.doi.org/10.3390/atmos11060594.

Castro, M.J., Ortega, S., de la Asunción, M., Mantas, J.M., Gallardo, J.M., 2011. GPU computing for shallow water flow simulation based on finite volume schemes. Comptes Rendus Mécanique 339 (2–3), 165–184.

Castro-Orgaz, O., Hager, W.H., 2017. Vertically Integrated Non-hydrostatic Free Surface Flow Equations. Springer International Publishing, Cham, pp. 17–79. http://dx.doi.org/10.1007/978-3-319-47971-2_2.

Caviedes-Voullième, D., Morales-Hernández, M., Norman, M.R., Özgen-Xian, I., v, Serghei(serghei-swe), 2023. 1.0: a performance-portable high-performance parallel-computing shallow-water solver for hydrology and environmental hydraulics. Geosci. Model. Dev. 16 (3), 977–1008. http://dx.doi.org/10.5194/gmd-16-977-2023.

Chen, J., Hill, A.A., Urbano, L.D., 2009. A gis-based model for urban flood inundation. J. Hydrol. 373 (1–2), 184–192.

Chevrel, M.O., Labroquère, J., Harris, A.J., Rowland, S.K., 2018. Pyflowgo: An open-source platform for simulation of channelized lava thermo-rheological properties. Comput. Geosci. 111, 167–180. http://dx.doi.org/10.1016/j.cageo.2017.11.009.

Cunge, J.A., Holly, F.M., Verwey, A., 1989. Practical Aspects of Computational River Hydraulics. Pitman Pub. Inc.

Delestre, O., Lucas, C., Ksinant, P.-A., Darboux, F., Laguerre, C., Vo, T.-N.-T., James, F., Cordier, S., 2013. Swashes: a compilation of shallow water analytic solutions for hydraulic and environmental studies. Internat. J. Numer. Methods Fluids 72 (3), 269–300. http://dx.doi.org/10.1002/fld.3741.

Echeverribar, I., Brufau, P., García-Navarro, P., 2023a. Extension of a roe-type riemann solver scheme to model non-hydrostatic pressure shallow flows. Appl. Math. Comput. 440, 127642.

Echeverribar, I., Brufau, P., García-Navarro, P., 2023b. A fully eulerian two-layer model for the simulation of oil spills spreading over coastal flows. Phys. Fluids 35 (11), 116602. http://dx.doi.org/10.1063/5.0169493.

Escalante, C., Dumbser, M., Castro, M.J., 2019. An efficient hyperbolic relaxation system for dispersive non-hydrostatic water waves and its solution with high order discontinuous galerkin schemes. J. Comput. Phys. 394, 385–416.

Escalante, C., Morales de Luna, T., Castro, M.J., 2018. Non-hydrostatic pressure shallow flows: Gpu implementation using finite volume and finite difference scheme. Appl. Math. Comput. 338, 631–659.

Ferrari, Stefania, Saleri, Fausto, 2004. A new two-dimensional shallow water model including pressure effects and slow varying bottom topography. ESAIM: M2AN 38 (2), 211–234. http://dx.doi.org/10.1051/m2an:2004010.

Finaud-Guyot, P., Hétier, M., Rousseau, A., 2023. Numerical modeling of large debris transport during floods. working paper or preprint.

Gairola, A., Bitsuamlak, G., 2019. Numerical tornado modeling for common interpretation of experimental simulators. J. Wind Eng. Ind. Aerodyn. 186, 32–48. http://dx.doi.org/10.1016/j.jweia.2018.12.013.

García-Martínez, R., Flores-Tovar, H., 1999. Computer modeling of oil spill trajectories with a high accuracy method. Spill Sci. Technol. Bull. 5 (5), 323–330. http://dx.doi.org/10.1016/S1353-2561(99)00077-8.

Alón-Rojas, I., Wang, X., Fredj, E., 2019. A 3d numerical model to track marine plastic debris (TrackMPD): Sensitivity of microplastic trajectories and fates to particle dynamical properties and physical processes. Marine Poll. Bull. 141, 256–272. http://dx.doi.org/10.1016/j.marpolbul.2019.02.052.

Knijff, J.M.V.D., Younis, J., Roo, A.P.J.D., 2010. LISFLOOD: A GIS-based distributed model for river basin scale water balance and flood simulation. Int. J. Geogr. Inf. Sci. 24 (2), 189–212.

Martínez-Aranda, S., Murillo, J., García-Navarro, P., 2020. A robust two-dimensional model for highly sediment-laden unsteady flows of variable density over movable beds. J. Hydroinformatics 22 (5), 1138–1160. http://dx.doi.org/10.2166/hydro.2020.027.

Merritt, D.M., Wohl, E.E., 2002. Processes governing hydrochory along rivers: Hydraulics, hydrology, and dispersal phenology. Ecol. Appl. 12 (4), 1071–1087. http://dx.doi.org/10.1890/1051-0761(2002)012[1071:PGHAR]2.0.CO;2.

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.

Morales-Hernández, M., Lacasta, A., Murillo, J., Brufau, P., García-Navarro, P., 2015. A Riemann coupled edge (RCE) 1D-2D finite volume inundation and solute transport model. Environ. Earth Sci. 74 (11), 7319–7335.

Morgan, R., Quinton, J., Smith, R., Govers, G., Poesen, J., Auerswald, K., Chisci, G., Torri, D., Styczen, M., Folly, A., 1998. The European Soil Erosion Model (EUROSEM): Documentation and User Guide. Cranfield University.

Nordam, T., Kristiansen, R., Nepstad, R., van Sebille, E., Booth, A.M., 2023. A comparison of eulerian and lagrangian methods for vertical particle transport in the water column. In: Geoscientific Model Development. vol. 16, (18), pp. 5339–5363. http://dx.doi.org/10.5194/gmd-16-5339-2023.

Ooyama, K., 1969. Numerical simulation of the life cycle of tropical cyclones. J. Atmos. Sci. 26 (1), 3–40. http://dx.doi.org/10.1175/1520-0469(1969)026<0003: NSOTLC>2.0.CO;2.
Ortega-Moya, J., Martínez-Aranda, S., Fernández-Pato, J., García-Navarro, P., 2024. A vertically non-uniform temperature approach for the friction term computation in depth-averaged viscoplastic lava flows. J. Comput. Phys. 519, 113378. http://dx.doi.org/10.1016/j.jcp.2024.113378.

Peeters, F., Hofmann, H., 2015. Length-scale dependence of horizontal dispersion in the surface water of lakes. Limnol. Oceanogr. 60, 1917–1934.

Persi, E., Petaccia, G., Sibilla, S., 2018. Large wood transport modelling by a coupled eulerian–lagrangian approach. Nat. Hazards 91 (1), 59–74, iD: Persi2018.

Petaccia, G., Leporati, F., Torti, E., 2016. Openmp and cuda simulations of sella zerbino dam break on unstructured grids. Comput. Geosci. 20 (5), 1123–1132.

Pilechi, A., Mohammadian, A., Murphy, E., 2022. A numerical framework for modeling fate and transport of microplastics in inland and coastal waters. Marine Poll. Bull. 184, 114119. http://dx.doi.org/10.1016/j.marpolbul.2022.114119.

Poland, M.P., Anderson, K.R., 2020. Partly cloudy with a chance of lava flows: Forecasting volcanic eruptions in the twenty-first century. J. Geophys. Res.: Solid Earth 125 (1), e2018JB016974. http://dx.doi.org/10.1029/2018JB016974.

Rutherford, J., 1994. River Mixing. Wiley, Chichester.

Sætra, M.L., Brodkorb, A.R., 2012. Shallow water simulations on multiple gpus. In: Jónasson, K. (Ed.), Applied Parallel and Scientific Computing. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 56–66.

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.

Sun, X., Zhao, C., Zhang, Y., Chen, F., Zhang, S., Zhang, K., 2021. Physical model test and numerical simulation on the failure mechanism of the roadway in layered soft rocks. Int. J. Min. Sci. Technol. 31 (2), 291–302. http://dx.doi.org/10.1016/j.ijmst.2021.01.003.

Teng, J., Jakeman, A., Vaze, J., Croke, B., Dutta, D., Kim, S., 2017. Flood inundation modelling: A review of methods, recent advances and uncertainty analysis. Environ. Model. Softw. 90, 201–216.

Thielen, J., Bartholmes, J., Ramos, M.-H., de Roo, A., 2009. The european flood alert system. PArt 1: Concept Dev. Hydrol. Earth Syst. Sci. 13 (2), 125–140.

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.

Vallés, P., Morales-Hernández, M., García-Navarro, P., Roeber, V., Caviedes-Voullième, D., 2023. Enhancing flood analysis with a lagrangian transport modeling and serghei. In: Advances in Hydroinformatics—SimHydro 2023, vol. 1, Springer, pp. 387–405. http://dx.doi.org/10.1007/978-981-97-4072-7.

Vallés, P., Morales-Hernández, M., Roeber, V., García-Navarro, P., Caviedes-Voullième, D., 2025. Serghei v2.1: a lagrangian model for passive particle transport using a 2d shallow water model (serghei-lpt). EGUsphere 2025, 1–27. http://dx.doi.org/10.5194/egusphere-2025-722.

Xia, X., Liang, Q., Ming, X., 2019. A full-scale fluvial flood modelling framework based on a high-performance integrated hydrodynamic modelling system (hipims). Adv. Water Resour. 132, 103392. http://dx.doi.org/10.1016/j.advwatres.2019.103392.

Yamazaki, Y., Kowalik, Z., Fai Cheung, K., Depth-integrated, 2009. Non-hydrostatic model for wave breaking and run-up. Int. J. Numer. Methods Fluids 61, 473–497.

You, H., Tinoco, R.O., 2023. Turbulent coherent flow structures to predict the behavior of particles with low to intermediate stokes number between submerged obstacles in streams. Water Resour. Res. 59 (2), e2022WR032439. http://dx.doi.org/10.1029/2022WR032439.