# SERGHEI v2.1: a Lagrangian model for passive particle transport using a two-dimensional shallow water model (SERGHEI-LPT)

Pablo Vallés, Mario Morales-Hernández, Volker Roeber, Pilar García-Navarro, Daniel Caviedes-Voullième

► To cite this version:

Pablo Vallés, Mario Morales-Hernández, Volker Roeber, Pilar García-Navarro, Daniel Caviedes-Voullième. SERGHEI v2.1: a Lagrangian model for passive particle transport using a two-dimensional shallow water model (SERGHEI-LPT). Geoscientific Model Development, 2025, 10.5194/gmd-18-7399-2025. doi:10.53196/gmd-18-7399

HAL Id: hal-05319676

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

Submitted on 17 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>

# SERGHEI v2.1: a Lagrangian model for passive particle transport using a two-dimensional shallow water model (SERGHEI-LPT)

Pablo Vallés¹,², Mario Morales-Hernández¹, Volker Roeber², Pilar García-Navarro¹, and Daniel Caviedes-Voullième³,⁴

¹I3A, University of Zaragoza, Zaragoza, Spain

²E2S Chair HPC-Waves, University of Pau, Anglet, France

³Institute of Bio- and Geosciences: Agrosphere (IBG-3), Forschungszentrum Jülich, Jülich, Germany

⁴Simulation and Data Lab. Terrestrial Systems, Jülich Supercomputing Centre (JSC), Forschungszentrum Jülich, Jülich, Germany

**Correspondence:** Daniel Caviedes-Voullième (d.caviedes.voullieme@fz-juelich.de)

Received: 15 February 2025 – Discussion started: 15 April 2025

Revised: 19 June 2025 – Accepted: 22 August 2025 – Published: 17 October 2025

**Abstract.** This paper presents a Lagrangian model for particle transport driven by a two-dimensional (2D) shallow water model, assuming that the particles have negligible mass and volume, are located at the free surface, and without interactions between them. Particle motion is based on advection and turbulent diffusion, which is added using a random-walk model. The equations for particle advective transport are solved using the flow velocity provided by a 2D shallow water solver and an online first-order Euler method, an online fourth-order Runge-Kutta method and an offline fourth-order Runge-Kutta method. The primary objective of this work is to present the capabilities of the new Lagrangian particle transport (LPT) model, while also providing an analysis of the accuracy and computational efficiency of the numerical schemes and their implementation for particle transport. To verify the accuracy and computational cost, several test cases inspired by laboratory set-ups are simulated. In this analysis, the Euler online method provides the best compromise between accuracy and computational efficiency. Finally, a localized precipitation event in the Arnás catchment is simulated to test the model’s capability to represent particle transport in overland flow over irregular topography.

## 1 Introduction

Over the past decade, floods have emerged as the most destructive natural disaster worldwide (Wallemacq et al., 2015; CRED and UNDRR, 2020), resulting in substantial and progressively increasing economic losses (Formetta and Feyen, 2019; Ripple et al., 2020). Projections suggest that in the coming decades, both economic losses and the number of individuals affected by floods are expected to rise (Dottori et al., 2018; OECD, 2016). This situation forces institutions to develop strategies and tools to predict and mitigate the damage caused by floods. Potential solutions range from administrative policies based on practical actions (United Nations, 2015; Olcina et al., 2016; Martin et al., 2010) to the development of computational prediction tools (Bates et al., 2023; Thielen et al., 2009; Chen et al., 2009). In recent years, these predictive tools, which rely on numerical simulations, have seen significant advancements in both accuracy and computational efficiency (Sampson et al., 2015; Knijff et al., 2010; Lacasta et al., 2014). Some models are based on the two-dimensional (2D) shallow water equations (2D SWE) model, offering accurate results with high computational efficiency (Xia et al., 2019; Morales-Hernández et al., 2021). However, the majority of models used for flood forecasting in rivers do not account for the transport of objects, such as vehicles, waste containers and other urban flood drifters (Bayón et al., 2024), or wood, carried by the flow during flooding events. This omission reduces the complete-
ness of these tools, as object transport is a critical compo-
nent of floods and can cause obstructions in hydraulic struc-
tures within channels and streams, thereby exacerbating the
damage during such events (Bayón et al., 2024; Valero et al.,
2024; Lofty et al., 2024). Several studies have examined de-
bris accumulation near weirs, culverts and dams, highlight-
ing their significant impact on flood risk (Thomas and Nisbet,
2012; De Cicco et al., 2015). Moreover, the transport of ob-
jects by geophysical flows is related to other issues, such as
environmental accidents due to pollutant spills in rivers and
lakes (van Emmerik et al., 2022; Ivshina et al., 2015; Mellink
et al., 2024), and environmental impact studies concerning
seed dispersion by precipitation in catchments (Zamora and
Montagnini, 2007; Sánchez-Salas et al., 2017), among oth-
ers.

The transport of objects by fluid flows has been exten-
sively studied in recent years (Braudrick and Grant, 2000;
Xia et al., 2011; Martínez-Gomariz et al., 2020). These stud-
ies have led to the development of computational models for
simulating the transport of particles with negligible mass or
volume (Finaud-Guyot et al., 2023), as well as contaminants
and small objects such as microplastics (García-Martínez
and Flores-Tovar, 1999; Jalón-Rojas et al., 2019), sediments,
bedload transport (Zhao et al., 2024; Baharvand et al., 2023)
and macroscopic objects (Persi et al., 2018a, b). Debris trans-
port can be analyzed using either an Eulerian or a Lagrangian
framework (Nordam et al., 2023). The Lagrangian approach
primarily provides information on the pathways linking the
origin to the destination of individual particles. It may also
capture specific processes affecting debris, such as deposi-
tion, fragmentation and degradation, provided that these pro-
cesses are explicitly implemented. On the other hand, the
Eulerian method provides a different approach by treating
debris as a concentration, offering greater computational ef-
ficiency than the Lagrangian method when the number of
bodies is very high. However, the Eulerian description is
inherently dispersive and does not account for the specific
transport trajectories of particles, thereby overlooking certain
small-scale details (Cai et al., 2023).

In recent decades, numerous computational models have
been developed to simulate particle transport. However, the
majority of these models are designed for ocean and some-
times coastal scenarios with fixed dry–wet boundaries (Le-
breton et al., 2012; Liubartseva et al., 2018). Furthermore,
some models update particle positions only at specific time
intervals rather than at every time step, to reduce the high
computational cost (Finaud-Guyot et al., 2023). These so-
called offline methods are often used in coastal environments
where the flow evolves on time scales much longer than
the particle time step, making temporal interpolation feasi-
ble (Cucco et al., 2009; Fajardo-Urbina et al., 2023, 2024). It
is important to note that flooding and drying occur in many
coastal systems, such as estuaries or tidal flats, and offline
methods have still been successfully applied in such con-
texts. The inaccuracy introduced by not updating particle

positions at every step is often mitigated by using higher-
order schemes, such as a fourth-order Runge-Kutta method
(García-Martínez and Flores-Tovar, 1999). An alternative ap-
proach is to determine the temporal evolution of material by
solving the Eulerian advection-dispersion equation, thereby
obtaining the changes in debris concentration over time (Ba-
harvand et al., 2023; Schreyers et al., 2024; Portillo De Arbe-
loa and Marzadri, 2024). While these models effectively cap-
ture the behavior of large quantities of small debris, they are
not suitable for large objects or scenarios involving a small
number of material particles. In addition, some Lagrangian
models are implemented using three-dimensional hydrody-
namic frameworks (Jalón-Rojas et al., 2019; Pilechi et al.,
2022), which provide more detailed flow information com-
pared to 2D and one-dimensional models, albeit at the cost
of increased computational time.

In addition to advection, turbulent diffusion can play an
important role in the material transport, not only in water
(Merritt and Wohl, 2002; Molazadeh et al., 2024; Yang and
Foroutan, 2023), but also in other environments such as wind
(Horn et al., 2012). Thus, it is interesting to add to the La-
grangian model of material transport some kind of minimal
turbulence model to achieve greater realism and accuracy in
the numerical results (Jalón-Rojas et al., 2019; Liubartseva
et al., 2018).

Finally, to be effective, the computational model must be
both accurate and computationally efficient. This balance is
especially important in applications such as flood forecast-
ing, real-time decision support, environmental impact assess-
ments and large-scale scenario simulations, where timely and
reliable results are essential. In this context, depth-averaged
2D hydrodynamic models have proven to offer a favorable
trade-off, providing sufficient accuracy for many surface wa-
ter flow scenarios while incurring significantly lower com-
putational costs compared to fully three-dimensional models
(Vacondio et al., 2016; Echeverribar et al., 2019).

This work represents an initial step towards a Lagrangian
model for material transport in shallow water flows. The ma-
terials simulated in this study are particles with negligible
mass and volume that do not interact with each other, thus
representing a simplified case. The primary objective is to
present the capabilities of the new Lagrangian particle trans-
port (LPT) model, together with an accompanying analysis
of the accuracy and computational efficiency of the numeri-
cal scheme and its implementation for particle transport. The
accuracy and efficiency of an explicit Euler method are com-
pared with those of a fourth-order Runge-Kutta method. In
addition, these methods are evaluated both when applied at
every time step or iteration, and at specific iteration inter-
vals. Furthermore, particle dispersion is incorporated into
the LPT model equations using a random-walk model be-
cause the hydrodynamic model does not resolve turbulence.
All analyses are based on simulations of both analytical and
laboratory test cases and realistic scenarios. The model is
driven by an Eulerian hydrodynamic framework based on
the SWE to describe flow evolution. The implementation of these techniques is part of an ongoing initiative to create a modeling framework grounded in physically based hydrodynamics: the Simulation Environment for Geomorphology, Hydrodynamics, and Eohydrology in Integrated form (SERGHEI) model (Caviedes-Voullième et al., 2023). The code developed for this study is open-source and available at https://gitlab.com/serghei-model (last access: 15 February 2025).

The structure of the paper is as follows: first, the governing equations for the hydrodynamic model and for the LPT model are presented; second, the numerical schemes and their details are presented; then, analytical, test and realistic cases are simulated; finally, the conclusions of the present work are shown.

## 2 Governing equations

The flow equations are presented and subsequently, the equations corresponding to passive particle transport are discussed.

### 2.1 2D SWE

The hydrodynamic model characterizes surface flow through the hyperbolic 2D SWE system, which is based on mass and momentum conservation (Cunge et al., 1989):

$$ \frac{\partial h}{\partial t} + \nabla(hv) = r - i + e + S_w, \quad (1) $$

$$ \begin{aligned} \frac{\partial hu}{\partial t} + \frac{\partial}{\partial x} \left( u^2 + \frac{1}{2} g h^2 \right) + \frac{\partial(huv)}{\partial y} & = -ghS_{fx} - gh \frac{\partial z_b}{\partial x}, \end{aligned} \quad (2) $$

$$ \begin{aligned} \frac{\partial hv}{\partial t} + \frac{\partial}{\partial y} \left( v^2 + \frac{1}{2} g h^2 \right) + \frac{\partial(huv)}{\partial x} & = -ghS_{fy} - gh \frac{\partial z_b}{\partial y}, \end{aligned} \quad (3) $$

where $h$ is the water depth [L], $v = (u, v)$ is the velocity flow vector [L T⁻¹], $r$ is the precipitation rate [L T⁻¹], $i$ is the infiltration rate [L T⁻¹], $e$ is the evaporation rate [L T⁻¹] and $S_w$ denotes the flow rate [L T⁻¹] provided by sources or sinks. The gravitational acceleration is denoted by $g$ [L T⁻²], $z_b$ is the bed elevation [L], and $S_{fx}$ and $S_{fy}$ are the friction slope terms, defined as

$$ 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 (4) $$

where $n$ is Manning's roughness coefficient [L T⁻¹/³] (Arlington and Schneider, 1984).

### 2.2 Particle transport equations

In a Lagrangian model of material transport, material elements are entrained and transported based on hydrodynamic forces, computed with appropriate coefficients or following a kinematic approach. In this study, the transported particles act as passive tracers: they have negligible volume and mass, do not interact with each other and do not affect the flow field.

In the absence of turbulence, particle transport is driven by advection, and the particle trajectories are governed by the following system of equations that defines the position of each particle $\mathbf{x}_p = (x_p, y_p, z_p)$:

$$ \left\{ \begin{aligned} \frac{dx_p}{dt} &= u(\mathbf{x}_p), \\ \frac{dy_p}{dt} &= v(\mathbf{x}_p), \\ z_p &= h(\mathbf{x}_p) + z_b(\mathbf{x}_p). \end{aligned} \right. \quad (5) $$

As observed, the vertical position is unaffected by the flow velocity due to the depth-averaged SWE approximation, which neglects the vertical velocity component. Consequently, the particle is assumed to reside at the free surface, computed as the sum of the water depth $h$ and the bed elevation $z_b$. This assignment is not derived from the governing equations but serves primarily to maintain numerical robustness and for visualization purposes. Since the hydrodynamic model is vertically averaged, and thus does not resolve vertical flow structure, particles transported by it do not possess a true vertical coordinate in the physical sense. However, assigning them a position at the free surface ensures consistency with the surface flow and avoids issues arising from irregular bathymetry. Notably, this choice prevents numerical artifacts, such as particles unrealistically crossing obstacles or walls due to inconsistent vertical velocities. Moreover, because the advection velocity is evaluated at the horizontal location of each particle, aligning all particles to the free surface provides a coherent reference for computing motion in the horizontal plane.

Turbulence effects in a flow can be significant and, by modifying the velocity field, also influence particle transport (Merritt and Wohl, 2002). Therefore, if turbulence is not incorporated into the flow velocity, as it is in the hydrodynamic model of the SERGHEI framework (Caviedes-Voullième et al., 2023), it can be included in the particle velocity, as follows:

$$ \left\{ \begin{aligned} \frac{dx_p}{dt} &= u(\mathbf{x}_p) + u_{\text{disp}}(\mathbf{x}_p), \\ \frac{dy_p}{dt} &= v(\mathbf{x}_p) + v_{\text{disp}}(\mathbf{x}_p), \\ z_p &= h(\mathbf{x}_p) + z_b(\mathbf{x}_p). \end{aligned} \right. \quad (6) $$

where $v_{\text{disp}} = (u_{\text{disp}}, v_{\text{disp}}, 0)$ is the velocity induced by the dispersion, where it is considered a null dispersion velocity in the vertical coordinate because of the lack of the advection velocity in the vertical. In this study, a random-walk model is employed to simulate dispersion because of its simplicity
and its great balance between computational efficiency and accuracy (Jalón-Rojas et al., 2019). This model proposes the $v_{disp}$ as a function of the diffusivity coefficients $K_{hx}$ and $K_{hy}$ [$L^2 T^{-1}$] 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. Therefore, the standard anisotropic diffusion model (Morales-Hernández et al., 2019) is utilized to determine the diffusivity variables, with the following expressions:

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

where $\epsilon_L$, and $\epsilon_T$ are the longitudinal and transversal dispersion coefficients [$L^2 T^{-1}$] (Rutherford, 1994), respectively; and $|u^*|$ is the friction velocity:

$$|u^*| = n \sqrt{g \frac{v_x^2 + v_y^2}{(h(x_p, y_p))^{1/3}}}. \quad (8)$$

### 3 Numerical schemes

Equations (1), (2) and (3) are solved using an explicit upwind finite-volume scheme, which is based on the augmented Roe–Riemann solver. Further details about the numerical scheme can be found in (Echeverribar et al., 2019; Murillo and García-Navarro, 2010; Morales-Hernández et al., 2013), together with the high-performance computing (HPC) implementation within SERGHEI (Caviedes-Voullière et al., 2023). In this paper, several numerical schemes are explored to solve the particle equations systems.

#### 3.1 Numerical scheme for the LPT model

To obtain the temporal and spatial evolution of the particles, it is necessary to discretize the system of Eq. (5) when turbulence is not considered, or Eq. (6) when it is included. In passive particle transport, particles are considered to have negligible mass and volume, ensuring no interaction with each other or the flow, thereby maintaining unaltered flow properties such as density or viscosity. The choice of discretization method depends on the required accuracy. In this study, two methods are utilized to discretize the equations governing particle trajectories: the forward Euler method and the Runge–Kutta method.

##### 3.1.1 Explicit forward Euler method

The temporal evolution of particle positions is obtained by considering the particles as mathematical points denoted by $p$, which are advected in space $x$ following a velocity field $v$. In this model, dispersion is not calculated within the hydrodynamic module, potentially leading to unrealistic particle

transport if this phenomenon is not considered (Merritt and Wohl, 2002). Although the flow has no turbulence model, the dispersive effects of turbulence are modeled on the transport of particles. These effects are added using a random-walk model, which can generate sufficient dispersion to introduce turbulent motion (Rutherford, 1994). Including the turbulence terms by using a random-walk model, the discretization of the system of Eq. (6) is as follows:

$$ 
\left\{
\begin{aligned}
x_p^{n+1} &= x_p^n + \Delta t_p^n u_i^n + \Delta t_p^n u_{disp,p}^n \\
&= x_p^n + \Delta t_p^n u_i^n + R_{p,x}^n (2\sigma^{-1} K_{hx,i} n \Delta t_p^n)^{1/2}, \\
y_p^{n+1} &= y_p^n + \Delta t_p^n v_i^n + \Delta t_p^n v_{disp,p}^n \\
&= y_p^n + \Delta t_p^n v_i^n + R_{p,y}^n (2\sigma^{-1} K_{hy,i} n \Delta t_p^n)^{1/2}, \\
z_p^{n+1} &= h_i^n + z_{b,i},
\end{aligned}
\right.
\quad (9) $$

with the dispersion velocities defined by

$$ u_{disp,p}^n = \sqrt{\frac{2K_{hx,i}^n}{\sigma \Delta t_p^n}}, \quad v_{disp,p}^n = \sqrt{\frac{2K_{hy,i}^n}{\sigma \Delta t_p^n}}, \quad (10) $$

where the Eulerian velocity field $v_i^n$ is the 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$ and $R_{p,y}^n$ are random numbers, which follow a uniform distribution with mean zero and standard deviation $\sigma = 1$; $K_{hx,i}^n$ and $K_{hy,i}^n$ are the horizontal diffusivity in $\hat{x}$ and $\hat{y}$ components, respectively, as defined by Eq. (7).

#### 3.1.2 Runge–Kutta method

The primary limitation of the explicit Euler method is its first-order accuracy, which may result in significant deviations from reality unless the time step for the Lagrangian model is sufficiently small (Bennett and Clites, 1987). Therefore, a higher-order approximation is necessary to capture the spatial variation of the velocity field when a large time step is computed by the numerical scheme. For this reason, a fourth-order Runge–Kutta method is presented for the advection formulation:

$$ 
\begin{aligned}
\mathbf{x}_p^{n+1} = \mathbf{x}_p^n + \frac{\Delta t_p^n}{6} \left( \mathbf{v}(\mathbf{x}_p, t^n) + 2\mathbf{v} \left( \mathbf{x}_p + \frac{1}{2}a, t^{n+1/2} \right) \right. \\
\left. + 2\mathbf{v} \left( \mathbf{x}_p + \frac{1}{2}b, t^{n+1/2} \right) + \mathbf{v} \left( \mathbf{x}_p + c, t^{n+1} \right) \right),
\end{aligned}
\quad (11) $$

where

$$ 
\begin{aligned}
a &= \Delta t_p^n \mathbf{v}(\mathbf{x}_p, t^n), & b &= \Delta t_p^n \mathbf{v}(\mathbf{x}_p + \frac{1}{2}a, t^{n+1/2}), \\
c &= \Delta t_p^n \mathbf{v}(\mathbf{x}_p + \frac{1}{2}b, t^{n+1/2}).
\end{aligned}
\quad (12) $$
<figure><img src="image_3.png" /><figcaption>Figure 1. The two possible strategies for the trajectory of particles in absence of obstacles (a) and the same strategies when there is an obstacle (gray cell) (b).</figcaption></figure>

This method requires the velocity at intermediate times ($t^{n+1/2}$), which is not directly calculated by the numerical scheme. Therefore, a fourth-order interpolation is used to obtain it:

$$v^{n+1/2}(\mathbf{x}_p) = \frac{5}{16}v^{n+1}(\mathbf{x}_p) + \frac{15}{16}v^n(\mathbf{x}_p) - \frac{5}{16}v^{n-1}(\mathbf{x}_p) + \frac{1}{16}v^{n-2}(\mathbf{x}_p). \quad (13)$$

Comparing the equations of both methods, it is evident that the higher accuracy of the Runge-Kutta method entails a higher computational cost per time step than the explicit Euler method, as the Runge-Kutta method requires the velocity vector at four different times and the calculation of additional terms to determine the trajectory evolution.

## 3.2 Trajectory algorithm

The transport of the particle is governed by the advection. The trajectory can be determined directly from the initial cell to the final cell (García-Martínez and Flores-Tovar, 1999; Bennett and Clites, 1987) (see red dashed trajectory in Fig. 1), or the particle can be moved cell by cell (see green dotted trajectory in Fig. 1). The latter approach is compulsory for addressing potential dry-wet interface issues and for preventing particles from traversing obstacles, as shown in Fig. 1b.

Thus, the algorithm that loops over the particles must be carefully designed to ensure that each particle traverses each computational cell discretely, dynamically adopting the local velocity and updating its position based on the discrete flow velocities within each cell.

For that, let us consider a single particle $p$ defined by its position $\mathbf{x}_p$ in the same position of the point $p = (p_x, p_y)$, which is inside cell $i$, with some velocity $v_i$. Let $d$ be the distance between the particle and the point at which the particle would leave the cell when following the velocity within the cell, as shown in Fig. 2a. In this figure, $\mathbf{A} = (A_x, A_y)$ and $\mathbf{q} = (q_x, q_y)$ are the lower-left and lower-right corners (re-

<figure><img src="image_4.png" /><figcaption>Figure 2. (a) Representation of a point $p$ (particle position) in a generic cell and (b) a generic representation of the time steps for 2D model and particles.</figcaption></figure>

spectively) of the cell, $b$ is the distance between $q$ and the intersection edge point and $\hat{e}$ is the unit vector of the Cartesian $y$-axis. The distance $d$ is the minimum distance in the direction of the velocity vector $v$ from the point $p$ until it crosses an edge that separates one computational cell from the next. This distance is measured along the line defined by the unit vector $\hat{v}$.

To guarantee that the particle is not advected beyond cell $i$ with the velocity $v_i$, it follows that a geometric restriction for the time step $\Delta t_p$ of particle $p$ is given by

$$\Delta t_p \le \frac{d}{||v_i||}. \quad (14)$$

This is fundamentally a Courant-Friedrichs-Lewy condition for the particle advection. An additional restriction applies for $\Delta t_p$ related to the hydrodynamic time-step size $\Delta t$:

$$\Delta t_p \le \Delta t, \quad (15)$$

ensuring the particle model cannot advance further in time than the hydrodynamic field. When $\Delta t_p < \Delta t$, particles must evolve through the full hydrodynamic time step via an iterative sub-cycling process. This involves calculating successive particle time steps $\Delta t_{p,m}$ until the total advancement matches the hydrodynamic time step (Fig. 2b):

$$\Delta t = \sum_{m=1}^{M} \Delta t_{p,m}, \quad (16)$$

where $M$ is the number of particle sub-steps within the hydrodynamic time step). This approach can produce situations like those in Fig. 3, where a particle traverses one, two or (at most) three cells during a single hydrodynamic time step $\Delta t$. Since particles can travel through at most three cells, the number of subdivisions $M$ satisfies $1 \le M \le 3$. The value of $M$ is determined individually for each particle based on the flow properties in the cells it traverses: $M = 1$ if the particle remains within the initial cell (Fig. 3a), $M = 2$ if it crosses into a neighboring cell (Fig. 3b), and $M = 3$ if it travels through three cells (Fig. 3c), though this case is uncommon
as it requires a specific combination of particle location and velocity field characteristics.

The complete formulation for the particle time-step restrictions is

$$ \Delta t_{p, m} \le \min \left( \frac{d}{\|\mathbf{v}_i\|}, \Delta t - \sum_{m=1}^{M-1} \Delta t_{p, m} \right), \quad (17) $$

where i refers to the current cell containing the particle, which may change during the hydrodynamic time step $\Delta t$.

The time-step restriction problem is therefore reduced to determining the distance d to the edge through which the particle will exit the cell.

The value of d is calculated as follows:

1. Given a line with direction $\hat{v} = v/|v|$ passing through particle point p, find its intersection with the lines defining the edges, and compute the distance d from p to the four intersection points: north, south, east and west intersections (see Fig. 2):

$$ p + d\hat{v} = q + b\hat{e}. \quad (18) $$

Considering only the x-coordinate:

$$ p_x + d\hat{v}_x = q_x, \quad (19) $$

and therefore, if $\hat{v}_x$ is not null,

$$ d = \frac{q_x - p_x}{\hat{v}_x} = \frac{A_x + \Delta x - p_x}{\hat{v}_x}, \quad (20) $$

where $\Delta x$ is the cell length in the x-coordinate, which, being a square cell, is equivalent to the length in the y-coordinate. For cases where $\hat{v}_x$ is null or negative, the distance d calculated using Eq. (20) does not yield the correct intersection. For this reason, three additional distances d are calculated by employing a procedure analogous to that used for deriving Eq. (20):

$$ d_2 = \frac{A_x - p_x}{\hat{v}_x} \quad \text{if and only if} \quad \hat{v}_x \neq 0, \quad (21) $$

$$ d_3 = \frac{A_y + \Delta x - p_y}{\hat{v}_y} \quad \text{if and only if} \quad \hat{v}_y \neq 0, \quad (22) $$

$$ d_4 = \frac{A_y - p_y}{\hat{v}_y} \quad \text{if and only if} \quad \hat{v}_y \neq 0. \quad (23) $$

2. Discard negative distances as these are antiparallel to the velocity vector.

3. The minimum positive distance d corresponds to the correct intersection.

Having computed d, it is possible to evaluate the particle sub-cycling time step (Eq. 17), and with that time step advect the particle with Eq. (6). This algorithm is used in both the Euler and Runge-Kutta methods.

Thus, the temporal evolution of the x- and y-coordinates of particles are obtained by the proposed algorithm. With respect to the z-coordinate, assuming a passive particle transport with negligible mass and volume for the particles, the value of the water surface level of the cell where the particle is located is imposed, such that it is floating in it. This approach and the use of the algorithm automatically resolve issues such as: (i) the dry-wet problem, provided the numerical method accurately addresses these conditions for the flow; and (ii) wet-wet problems involving topographical changes, where a particle floating at an arbitrary vertical position within a cell may encounter an adjacent cell with a higher elevation, preventing its movement into that cell despite both cells containing water. This ensures that particle movement is physically feasible and adheres to the elevation constraints of the cells.

Introducing a turbulent dispersion term can cause issues in wet-dry cases, as a particle could find itself in a wet cell due to advection and then reach a dry cell when the dispersion term is applied. Therefore, a restriction is imposed: if the destination cell (determined by adding the dispersion term) is dry, the dispersion term is not applied. This restriction is implemented for several reasons:

- It is physically inconsistent for a particle with negligible mass to enter a dry cell. The underlying principle is that particles with significant mass could use inertia to traverse dry regions. Therefore, this restriction specifically applies to massless or negligible-mass particles.

- The dispersion term can relocate the particle horizontally across the surface, and if this displacement places it in a dry cell with a higher elevation than its previous location, the particle undergoes an unphysical “jump” in its trajectory. This arises because dispersion terms depend on friction velocity without properly accounting for flow direction constraints.

### 3.3 Temporal integration modes

Depending on the desired accuracy of the results, different time intervals are proposed to update the evolution of the particles. When the particle advection is computed on every time step of the hydrodynamic solver (from n to n + 1), this configuration is referred to as the online mode. In contrast, when particle positions are updated every n/100 hydrodynamic time steps (from n to n + n/100), the configuration corresponds to the offline mode. The online mode is expected to yield higher accuracy than the offline mode, but potentially with a higher computational cost due to more frequent advection updates. The explicit Euler method is unsuitable for offline mode due to its limited accuracy, and is therefore applied.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Three possible situations for a particle: traveling through 1 cell (a), 2 cells (b), or 3 cells (c).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Conceptual diagram for the SERGHEI modules and their implementation status.</figcaption></figure>

plied exclusively in the online configuration. Conversely, the Runge–Kutta method demonstrates sufficient accuracy to be employed in both online and offline modes.

## 3.4 Implementation in SERGHEI framework

The LPT model has been implemented in the SERGHEI framework, as one of the several modules that have been implemented in this framework (see Fig. 4). The SERGHEI framework is a modular and open-source simulation system designed to target performance portability across HPC platforms and standard workstations (Caviedes-Voulième et al., 2023). It focuses on future-proofing and sustainability while consolidating over a decade of mature numerical techniques for shallow water modeling. SERGHEI facilitates a modular, HPC-ready community model, enabling broader collaboration. The framework moves beyond classical fluvial flooding and hydrograph generation to address applications related to landscape function, transport and integration into Earth system models.

The LPT model has been developed along the same lines as the other modules in SERGHEI. In this way, this model has been implemented in C++ and using the Kokkos performance-portability layer, which allows the simulation of the model using any GPU architecture. As can be seen in Fig. 5, the LPT model is updated after updating the hydro-

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Schematic of the numerical calculation using the LPT module.</figcaption></figure>

dynamic model, as it needs the updated information of the conserved variables. The LPT model is currently being implemented to support distributed computations both on multiple CPU nodes (note that shared memory CPU parallelization is achieved via OpenMP) and on multi-GPU systems, following the approach used in other SERGHEI modules.

## 4 Results

The validation of the LPT model is conducted through the simulation of several test cases with increasing complexity, assessing both accuracy and computational efficiency. The
simulations are run using a GPU NVIDIA GeForce RTX 3060.

## 4.1 Steady circular vortex case

Initially, a steady circular vortex test case is simulated to validate the model's accuracy. The vortex has a diameter of 10 m and is centered in the middle of the domain. Its velocity field increases radially from zero at the center to a maximum of 100 m s⁻¹ at the outer edge. A steady-state configuration is assumed, meaning that particles released into the flow should maintain a constant radius throughout their trajectories (Finaud-Guyot et al., 2023; García-Martínez and Flores-Tovar, 1999). To ensure reproducibility, the particle is released at a distance of 10 m from the vortex center. The circular vortex is simulated with varying cell sizes, ensuring the particle completed four full revolutions for consistency in error measurement. The vortex radius is consistently set at R = 100 m, with a water depth of h = 0.1 m and no roughness (n = 0). For this test, the hydrodynamic solver is turned off, so that the Eulerian flow fields (water depth and velocity) remained constant and are not evolved over time, also resulting in a constant time step Δt throughout the simulation. Turbulent diffusion terms are excluded from the particle trajectory equations because this is an analytical case, and the objective is to guarantee correctness and robustness of the advection. The Runge-Kutta offline method is updated every five hydrodynamic time steps, a choice made to balance computational efficiency and accuracy. Preliminary tests with shorter update intervals (two to three time steps) yielded only marginal cost savings and limited improvement in accuracy, so the five-step update was selected to better illustrate this trade-off.

The errors for each grid resolution using different numerical configurations are shown in Fig. 6. The L₁-norm for the error at each time step is defined by

$$L_1(t^n) = |x_p^n - \hat{x}|, \qquad (24)$$

where $x_p^n$ is the particle position at the time $t^n$ and $\hat{x}$ is the analytical position. The L₁-norm integrated in time is the mean absolute error (MAE), which is obtained using

$$\text{MAE} = \frac{1}{N} \sum_{n=0}^{N} L_1(t^n), \qquad (25)$$

where N is the number of time measurements.

In addition, the L₂-norm at each time step is defined as

$$L_2(t^n) = (x_p^n - \hat{x})^2, \qquad (26)$$

while the root-mean-squared error (RMSE) is computed as

$$\text{RMSE} = \sqrt{\frac{1}{N} \sum_{n=0}^{N} L_2(t^n)}. \qquad (27)$$

Figure 6 shows that, under steady flow, the error of the Runge-Kutta method is similar to that of the Euler method. The algorithm implemented provides a reduced error for the Euler method because of the calculus of the trajectory in each cell (see Fig. 1) and in each time step. Moreover, the error of all methods depends strongly on the discretization of the domain, with lower errors for smaller cell sizes. Figure 6 also shows the error reduction in the same order as the Euler method as the cell size is reduced (order one). However, for the Runge-Kutta method (order four), the error is not observed to reduce with that order. This is mainly because it is a steady case and the value of the time step is fixed.

## 4.2 Transient test case

Following the accuracy assessment, an analysis of computational efficiency is conducted using a well-known test case (see Fig. 7a), previously simulated and compared against measurements using the SERGHEI framework (Caviedes-Voullième et al., 2023). The details of the experimental setup are described in (Soares-Frazão and Zech, 2008), where the domain dimensions are L = 36 m in length, B = 3.6 m in width. As the initial condition, an initial water depth of 0.4 m in the reservoir region and 0 m in the right-hand region are imposed, modifying the value of the initial condition in the right-hand region compared to the original experiment to check the accuracy of the model when wet-dry situations occur. The comparative analysis is based on increasing the number of particles to evaluate the model computational cost for each configuration. Moreover, this case is used as a validation of the model's ability to correctly resolve wet-dry situations. Reflective boundary conditions are applied everywhere, maintaining a constant number of particles throughout the simulation. The Runge-Kutta offline method was updated using different time periods (5Δt, 10Δt, 20Δt and 50Δt, where Δt is the hydrodynamic time step) to analyze the influence of this variable in the temporal evolution of the particle position.

The computational costs for the different configurations with different particle numbers are listed in Table 1, where only the Runge-Kutta offline computational costs updated every five time steps are shown because it is the less computationally efficient case for the simulated Runge-Kutta offline combinations. The computational time for this simulation without particles was 12.15 s, and an overhead ratio relative to the hydrodynamic-only run is computed for each case. The forward Euler method consistently demonstrated shorter computational times compared to other methods across all configurations. In addition, four states during the simulation are depicted in Fig. 7, illustrating realistic particle trajectories simulated with the Euler scheme without traversing dry cells or climbing over buildings. However, when the Runge-Kutta offline is updated using long time periods, the temporal evolution of particle position can be illogical, as shown in Fig. 8, where some particles can be seen climbing over
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Errors simulating the steady circular vortex case using different numerical methods, shown in logarithmic scales: MAE (a) and RMSE (b). A first-order convergence slope is included for reference.</figcaption></figure>

buildings. Finally, the differences between the methods' results are shown in Fig. 9, where the $L_2$-norm per particle is derived from the differences between models, taking the Runge–Kutta online method as the basis for comparison. The error between the Euler and Runge–Kutta methods is lower than the errors between the Runge–Kutta and Runge–Kutta offline methods, indicating that the Euler method provides results more similar to the Runge–Kutta method while maintaining higher computational efficiency. In addition, as the update time for the offline Runge–Kutta method increases, the difference between the online Runge–Kutta and the offline Runge–Kutta increases too, due to the fact that being a transient case, not updating frequently can lead to a high error.

## 4.3 Channel with cavities

A longitudinal channel with symmetrically shaped cavities was simulated, using a mesh with a resolution of 2 × 2 cm (see Fig. 10a). The objective is to check whether a symmetrical distribution of particles is achieved across the upper and lower cavities (Vallés et al., 2023); 10⁵ particles were simulated with an initial distribution as shown in Fig. 10a. Given the best compromise between accuracy and computational efficiency in the last test cases, the Euler method was only employed for this simulation. The temporal evolution of the particle position with and without turbulence was compared (see Fig. 10), with the evolution of particle numbers in each pair of upper and lower cavities depicted in Fig. 11. The presence of turbulence resulted in greater dispersion in particle distribution and more realistic trajectories compared to simulations without turbulence. This is further evidenced in Fig. 11, where the number of particles within a cavity is higher when turbulence is considered. The symmetry between the upper and lower cavities (see Fig. 11), regardless of whether turbulence is added or not, provides a further argument for the accuracy of the Lagrangian model.

An important observation is that with coarser grid resolutions, the velocity field may not be sufficiently resolved to drive particles into the cavities. To explore this, the same simulation was repeated with and without dispersion using a different grid resolution: 20 cm. As shown in Figs. 12 and 13, the coarser mesh failed to guide particles into several cavities, with some receiving no particles at all. In contrast, the finer resolution (see Fig. 11) captured the expected cavity entry more accurately. This highlights a key limitation of coarse Eulerian velocity fields, which can be partially mitigated by including the random-walk model: it enhances subgrid-scale transport and improves the physical plausibility of particle trajectories. Although this test does not constitute formal model verification, the observed symmetry in particle distributions – particularly under the influence of turbulence – provides soft validation that the Lagrangian model behaves in a consistent and physically reasonable way, without introducing spurious asymmetric artifacts.

## 4.4 Example application: surface runoff in the Arnás catchment

The implemented Lagrangian model may have applications in ecohydrology, such as seed, pollutant or nutrient transport. Consequently, a case involving a precipitation event over the Arnás catchment (see Fig. 14a) was simulated to show the potential application of this model. This catchment is located in the Borau valley (northern Spanish Pyrenees). In the last decades, this region has been an area of exhaustive study of different hydrological processes (Lana-Renault et al., 2007; García-Ruiz et al., 2005). Physical processes such as rainfall, infiltration, runoff generation and gully formation have been simulated in detail using various methods (López-Barrera et al., 2011; Fernández-Pato et al., 2016; Vallés et al., 2024). The primary objective of simulating this case was to assess its potential utility in environmental applications. Information derived from conservative tracers and field experiments measuring isotopes in water has been widely used in hydrology
<table><caption>**Table 1.** Computational cost for different number of particles using the numerical methods implemented and the increase ratio between simulation with particles and without particles, for the dam break test case.</caption><thead><tr><td>Numerical method</td><td>Number of particles</td><td>Computational cost (s)</td><td>Increase ratio</td></tr></thead><tbody><tr><td>Euler</td><td>10<sup>4</sup></td><td>13.15</td><td>1.08</td></tr><tr><td>Euler</td><td>10<sup>5</sup></td><td>21.22</td><td>1.75</td></tr><tr><td>Euler</td><td>10<sup>6</sup></td><td>66.59</td><td>5.48</td></tr><tr><td>Runge-Kutta</td><td>10<sup>4</sup></td><td>30.00</td><td>2.47</td></tr><tr><td>Runge-Kutta</td><td>10<sup>5</sup></td><td>121.93</td><td>10.03</td></tr><tr><td>Runge-Kutta</td><td>10<sup>6</sup></td><td>1043.56</td><td>85.89</td></tr><tr><td>Runge-Kutta offline 5Δt</td><td>10<sup>4</sup></td><td>18.89</td><td>1.55</td></tr><tr><td>Runge-Kutta offline 5Δt</td><td>10<sup>5</sup></td><td>60.93</td><td>5.01</td></tr><tr><td>Runge-Kutta offline 5Δt</td><td>10<sup>6</sup></td><td>474.24</td><td>39.03</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Initial condition of the transient test case (a). States at $t = 5$ s (b), at $t = 10$ s (c), at $t = 15$ s (d), and at $t = 20$ s (e) for the Euler method with $10^5$ particles in the dam break case.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** State with particles climbing over buildings for the Runge-Kutta offline simulation with time step $50Δt$ (left), and zoom of the building area (right).</figcaption></figure>

to characterize hydrological pathways and connectivity, with ideas such as transit or residence times, and travel time distributions (Benettin et al., 2022). Here, a minimal example is presented to demonstrate how analogous ideas can be investigated within the proposed computational modeling framework.

The catchment was discretized at 5 m resolution, for a total of $1.1 \times 10^5$ computational cells; $10^4$ particles were randomly placed in the domain as initial condition for the LPT model, and the catchment was initially dry for the shallow water solver. As in the previous test case, the Euler method was exclusively used for the simulation of the Arnás case, based on the optimal balance between accuracy and computational efficiency observed in the initial test cases. Two real precip-
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Temporal evolution of the $L_2$-norm per particle obtained from the difference between methods results.</figcaption></figure>

iteration cases (see Fig. 14b and c) were simulated to analyze the particle trajectory behavior under different hyetographs. In both events, particles are mostly transported through the main gullies. The run times for the simulations of events 1 and 2 were 317.75 and 177.02 s, respectively, corresponding to run-time increase ratios of 2.39× and 1.10× compared to the same simulations without particles.

In addition, the total travel distance per particle, the particle travel time (travel time), and the time during which the particle was at rest (stop time) are shown in Figs. 15 and 16 for each precipitation event. These figures provide a statistical overview of particle dynamics during the runoff events. Physically, these metrics offer insights into the hydrological processes at play. The travel distance (displayed in Fig. 14a) reflects the efficiency of the drainage network of the basin to transport the debris, which can be seen in Fig. 15a and b, where most particles travel a distance very similar to the initial one. Moreover, Fig. 15a shows that a subset of particles from event 2 travels much farther than the rest (around 10 km). This is due to the properties of this precipitation event, which features less intense but longer rainfall than event 1. This extended duration causes particles to travel for longer periods (with shorter stop times), as also observed in Fig. 15d. Furthermore, this histogram reveals regions where particles tend to accumulate or remain at rest for extended periods, typically due to reduced flow velocities or topographical barriers, which may indicate zones of storage or potential sediment deposition. However, as shown in Fig. 16, the set of particles with high travel distances does not influence the last hours of the temporal evolution of the total travel distance (see Fig. 16) because in both events the majority of the particles leave the domain. However, the difference between the events is higher at the beginning of the case, having a big difference between both hyetographs in the first hours. In event 1, the longer travel times observed for most particles

in Fig. 15c, where particles tend to remain approximately at rest for extended periods, suggest intermittent flow conditions, possibly due to intermittent rainfall, leading to temporary particle deposition. In contrast, event 2 exhibits more uniform movement and shorter stop times, probably due to less intermittent rainfall that fosters persistent particle transport. Thus, these histograms not only quantify particle behavior but also provide insights into the underlying hydrological and geomorphological processes shaping the catchment's response to precipitation events.

## 5 Conclusions

In this work, a Lagrangian model for particle transport implemented within the SERGHEI framework is presented. This model utilizes an algorithm that leverages information from each cell occupied by the particles and is solved using three numerical schemes: an online first-order Euler method, an online fourth-order Runge-Kutta method, and an offline fourth-order Runge-Kutta method.

The convergence and accuracy of these solvers were first assessed using a steady circular vortex benchmark case. All methods yielded comparable error metrics and demonstrated proper grid convergence behavior. In particular, the Euler and fourth-order Runge-Kutta online schemes produced nearly identical results under steady conditions, highlighting the accuracy of the Euler method despite its lower order.

Next, a transient test case without dispersive terms was examined to evaluate computational efficiency and the handling of dry-wet interfaces. Results showed that the Euler method was the most computationally efficient, with relatively small differences in accuracy compared to the fourth-order Runge-Kutta schemes. Furthermore, as shown in Fig. 9, the offline fourth-order Runge-Kutta method produced incorrect trajectories due to its inability to account for dry cells during interpolation. The high accuracy of the Euler scheme is attributed to the sufficiently small time steps used in the simulations. While fourth-order Runge-Kutta online may offer slightly improved accuracy in transient flows with fewer particles, its computational cost increases significantly with particle count; up to 17 times more than Euler in cases involving 10⁵ particles or more. Given that this model is designed for large-scale applications, where millions of particles may be simulated, the Euler method offers the best trade-off between accuracy and efficiency, making it the preferred scheme for practical use.

To investigate the effect of turbulence, a channel with multiple cavities, with and without turbulence, was simulated. The results illustrate that incorporating turbulence significantly alters particle trajectories, resulting in more realistic transport dynamics. However, as it is shown in previous work (Vallés et al., 2024), the addition of turbulence implies a relevant increase in the computational cost. In addition, a comparison with a coarser grid resolution demonstrated that in-
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Initial condition for the cavities test case with 10⁵ particles with the numeration of the polygons (a). States with zoom at *t* = 50 s without turbulence (b) and with turbulence (c), and at *t* = 100 s without turbulence (d) and with turbulence (e).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Temporal evolution of the number of particles inside of the different cavities (cavities 1 and 5 in a, cavities 2 and 6 in b, cavities 3 and 7 in c, and cavities 4 and 8 in d).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** States with zoom at *t* = 50 s without turbulence (a) and with turbulence (b), and at *t* = 100 s without turbulence (c) and with turbulence (d) for a resolution of 20 cm.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Temporal evolution of the number of particles inside of the different cavities (cavities 1 and 5 in a, cavities 2 and 6 in b, cavities 3 and 7 in c, and cavities 4 and 8 in d) for a resolution of 20 cm.</figcaption></figure>

sufficient spatial resolution in the Eulerian velocity field may prevent particles from entering cavity regions. This limitation can be partially mitigated by including the random-walk dispersion model, which enhances subgrid-scale transport and contributes to more physically consistent trajectories. The observed symmetry in particle distributions, especially under turbulent conditions, provides further support for the model’s qualitative reliability.

Finally, the model was applied to a realistic case study of the Arnás catchment. The output generated by the model offers significant insights for environmental applications. The results elucidate the mechanisms of particle transport during precipitation events and identifies the primary gullies that fa-

cilitate this transport. Furthermore, this information provides critical insight into the dynamics of particle movement, revealing distinct patterns in both movement and stop times. These patterns reflect the influence of varying precipitation intensities and can help on identifying sediment mobilization and deposition processes within the catchment. This analysis underscores the potential of such models to enhance our understanding of different multiscale complex phenomena that take place in a catchment.

In summary, the Euler method delivers a strong balance between computational efficiency and numerical accuracy, particularly in large-scale simulations where performance is critical. The LPT model developed within the SERGHEI
<figure><img src="image_16.png" /><figcaption>**Figure 14.** 3D representation of the Arnás catchment, including the travel distance and the initial distance (a), and hyetographs of event 1 (b) and event 2 (c) for the Arnás case.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Histograms of: (a) the travel distance, (b) the travel distance normalized by the initial distance between the particle and the boundary condition position for both precipitation events, (c) the travel time, and (d) the stop time.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Temporal evolution of the total travel distance for both events.</figcaption></figure>

framework effectively reproduces realistic particle trajectories and can account for complex interactions such as dry-wet interfaces and turbulence. Future work will focus on optimizing the model to further reduce computational costs and implementing multi-GPU simulations to leverage the capabilities of the SERGHEI hydrodynamic model. Moreover, several enhancements will be incorporated into the LPT module to increase the realism of particle trajectories. These improvements will enable representation of pollutant transport (e.g., microplastics) and biological dispersal (e.g., seeds). In addition, incorporating particle mass, volume and inertia will allow modeling of vertical movement, deposition processes and macroscopic objects such as logs or urban debris in floods. However, accurate simulation of these transport phenomena requires development and validation of specific physical processes that are currently beyond the scope of the present model.

**Code and data availability.** SERGHEI is available through GitLab, at https://gitlab.com/serghei-model/serghei (last access: 15 February 2025), under a 3-clause BSD license. SERGHEI v2.1.0 was tagged as the first release of the LPT module at the time of submission of this paper. Static code for v2.1.0 can be found at https://doi.org/10.5281/zenodo.14871005 (Vallés et al., 2025b). This release does not include the Runge–Kutta solver. Static code for the experimental solution with Runge–Kutta integrators can be found at https://doi.org/10.5281/zenodo.14870918 (Vallés et al., 2025a). A repository containing test cases is available https://gitlab.com/serghei-model/serghei-tests/lpt (last access: 15 February 2025).

**Author contributions.** PV: investigation, methodology, software, validation, writing – original draft preparation. DCV: supervision, conceptualization, methodology, software, formal analysis, data cu-

ration, writing – review and editing. MMH: supervision, funding acquisition, conceptualization, methodology, formal analysis, writing – review and editing. VR: supervision, funding acquisition, writing – review and editing. PGN: funding acquisition, writing – review and editing.

**Competing interests.** The contact author has declared that none of the authors has any competing interests.

**Disclaimer.** Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

**Acknowledgements.** Pablo Vallés is funded by the UPPA-UNIZAR Research Grant PI-PRD/2022-03 and was previously supported by the Erasmus+ KA103 “IBERUSC” (2021-1-ES01-KA130-HED-000004265). The authors also acknowledge the Government of Aragón (Fluid Dynamics Technologies, T32_23R), the Universidad de Zaragoza (project JIUZ2023-IA-04), and Forschungszentrum Jülich for institutional support.

**Financial support.** This research has been supported by the Universidad de Zaragoza (grants PI-PRD/2022-03 and JIUZ2023-IA-04), the Université de Pau et des Pays de l’Adour (grant PI-PRD/2022-03), the Erasmus+ program (grant KA103 IBERUSC 2021-1-ES01-KA130-HED-000004265), the project PID2022-137334NB-I00 funded by MCIN/AEI/10.13039/501100011033 and ERDF/EU, and the Government of Aragón (grant T32_23R, Fluid Dynamics Technologies). Institutional support was also provided by Forschungszentrum Jülich.

The article processing charges for this open-access publication were covered by the Forschungszentrum Jülich.

**Review statement.** This paper was edited by Deepak Subramani and reviewed by two anonymous referees.

## References

Arcement, G. and Schneider, V.: Guide for Selecting Manning's Roughness Coefficients for Natural Channels and Flood Plains, no. 2339 in U.S. Geological Survey, Water-supply paper, https://pubs.usgs.gov/publication/wsp2339 (last access: 9 October 2025), 1984.

Baharvand, S., Ahmari, H., and Taghvaei, P.: Developing a Lagrangian sediment transport model for open channel flows, International Journal of Sediment Research, 38, 153–165, https://doi.org/10.1016/j.ijsrc.2022.09.003, 2023.
Bates, P. D., Savage, J., Wing, O., Quinn, N., Sampson, C., Neal, J., and Smith, A.: A climate-conditioned catastrophe risk model for UK flooding, Nat. Hazards Earth Syst. Sci., 23, 891–908, https://doi.org/10.5194/nhess-23-891-2023, 2023.

Bayón, A., Valero, D., and Franca, M. J.: Urban flood drifters (UFD): Identification, classification and characterization, Journal of Flood Risk Management, e13002, https://doi.org/10.1111/jfr3.13002, 2024.

Bennett, J. R. and Clites, A. H.: Accuracy of trajectory calculation in a finite-difference circulation model, Journal of Computational Physics, 68, 272–282, https://doi.org/10.1016/0021-9991(87)90058-1, 1987.

Benettin, P., Rodriguez, N. B., Sprenger, M., Kim, M., Klaus, J., Harman, C. J., van der Velde, Y., Hrachowitz, M., Botter, G., McGuire, K. J., Kirchner, J. W., Rinaldo, A., and McDonnell, J. J.: Transit Time Estimation in Catchments: Recent Developments and Future Directions, Water Resources Research, 58, https://doi.org/10.1029/2022wr033096, 2022.

Braudrick, C. A. and Grant, G. E.: When do logs move in rivers?, Water Resources Research, 36, 571–583, https://doi.org/10.1029/1999WR900290, 2000.

Cai, C., Zhu, L., and Hong, B.: A review of methods for modeling microplastic transport in the marine environments, Marine Pollution Bulletin, 193, 115136, https://doi.org/10.1016/j.marpolbul.2023.115136, 2023.

Caviedes-Voullième, D., Morales-Hernández, M., Norman, M. R., and Özgen-Xian, I.: SERGHEI (SERGHEI-SWE) v1.0: a performance-portable high-performance parallel-computing shallow-water solver for hydrology and environmental hydraulics, Geosci. Model Dev., 16, 977–1008, https://doi.org/10.5194/gmd-16-977-2023, 2023.

Chen, J., Hill, A. A., and Urbano, L. D.: A GIS-based model for urban flood inundation, Journal of Hydrology, 373, 184–192, 2009.

CRED and UNDRR (Centre for Research on the Epidemiology of Disasters and United Nations Office for Disaster Risk Reduction): Human Cost of Disasters: An overview of the last 20 years (2000–2019), Centre for Research on the Epidemiology of Disasters (CRED), Brussels, https://www.undrr.org/publication/human-cost-disasters-overview-last-20-years-2000-2019 (last access: 9 October 2025), 2020.

Cucco, A., Umgiesser, G., Ferrarin, C., Perilli, A., Canu, D. M., and Solidoro, C.: Eulerian and lagrangian transport time scales of a tidal active coastal basin, Ecological Modelling, 220, 913–922, https://doi.org/10.1016/j.ecolmodel.2009.01.008, 2009.

Cunge, J. A., Holly, F. M., and Verwey, A.: Practical aspects of computational river hydraulics, Pitman Pub. Inc., https://cir.nii.ac.jp/crid/1971149384816685571 (last access: 9 October 2025), 1989.

De Cicco, P., Enio, P., and Solari, L.: Flume experiments on bridge clogging by woody debris: The effect of shape of piers, International Association for Hydro-Environment Engineering and Research (IAHR), https://doi.org/10.13140/RG.2.1.4630.6648, 2015.

Dottori, F., Szewczyk, W., Ciscar, J.-C., Zhao, F., Alfieri, L., Hirabayashi, Y., Bianchi, A., Mongelli, I., Frieler, K., Betts, R. A., and Feyen, L.: Increased human and economic losses from river flooding with anthropogenic warming, Nature Climate Change, 8, 781–786, 2018.

Echeverribar, I., Morales-Hernández, M., Brufau, P., and García-Navarro, P.: 2D numerical simulation of un-

steady flows for large scale floods prediction in real time, Advances in Water Resources, 134, 103444, https://doi.org/10.1016/j.adwatres.2019.103444, 2019.

Fajardo-Urbina, J. M., Arts, G., Gräwe, U., Clercx, H. J. H., Gerkema, T., and Duran-Matute, M.: Atmospherically Driven Seasonal and Interannual Variability in the Lagrangian Transport Time Scales of a Multiple-Inlet Coastal System, Journal of Geophysical Research: Oceans, 128, e2022JC019522, https://doi.org/10.1029/2022JC019522, 2023.

Fajardo-Urbina, J. M., Liu, Y., Georgievska, S., Gräwe, U., Clercx, H. J., Gerkema, T., and Duran-Matute, M.: Efficient deep learning surrogate method for predicting the transport of particle patches in coastal environments, Marine Pollution Bulletin, 209, 117251, https://doi.org/10.1016/j.marpolbul.2024.117251, 2024.

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

Finaud-Guyot, P., Hétier, M., and Rousseau, A.: Numerical modeling of large debris transport during floods. HAL-Inria, https://inria.hal.science/hal-04009769 (last access: 9 October 2025), 2023.

Formetta, G. and Feyen, L.: Empirical evidence of declining global vulnerability to climate-related hazards, Global Environmental Change, 57, 101920, https://doi.org/10.1016/j.gloenvcha.2019.05.004, 2019.

García-Martínez, R. and Flores-Tovar, H.: Computer Modeling of Oil Spill Trajectories With a High Accuracy Method, Spill Science & Technology Bulletin, 5, 323–330, https://doi.org/10.1016/S1353-2561(99)00077-8, 1999.

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

Horn, S., Raabe, A., Will, H., and Tackenberg, O.: TurbSeed – A model for wind dispersal of seeds in turbulent currents based on publicly available climate data, Ecological Modelling, 237-238, 1–10, https://doi.org/10.1016/j.ecolmodel.2012.04.009, 2012.

Ivshina, I., Kuyukina, M., Krivoruchko, A., Elkin, A., Makarov, S., Cunningham, C., Peshkur, T., Atlas, R., and Philp, J.: Oil spill problems and sustainable response strategies through new technologies, Environmental Science: Processes and Impacts, 17, 1201–1219, https://doi.org/10.1039/c5em00070j, 2015.

Jalón-Rojas, I., Wang, X., and Fredj, E.: A 3D Numerical Model to Track Marine Plastic Debris (TrackMPD): Sensitivity of Microplastic Trajectories and Fates to Particle Dynamical Properties and Physical Processes, Marine Pollution Bulletin, 141, 256–272, https://doi.org/10.1016/j.marpolbul.2019.02.052, 2019.

Knijff, J. M. V. D., Younis, J., and Roo, A. P. J. D.: LISFLOOD: A GIS-based distributed model for river basin scale water balance and flood simulation, International Journal of Geographical Information Science, 24, 189–212, 2010.

Lacasta, A., Morales-Hernández, M., Murillo, J., and García-Navarro, P.: An optimized GPU implementation of a 2D free surface simulation model on unstructured meshes, Advances in Engineering Software, 78, 1–15, 2014.
Lana-Renault, N., Latron, J., and Regués, D.: Streamflow response and water-table dynamics in a sub-Mediterranean research catchment (Central Pyrenees), Journal of Hydrology, 347, 497–507, https://doi.org/10.1016/j.jhydrol.2007.09.037, 2007.

Lebreton, L.-M., Greer, S., and Borrero, J.: Numerical modelling of floating debris in the world’s oceans, Marine Pollution Bulletin, 64, 653–661, https://doi.org/10.1016/j.marpolbul.2011.10.027, 2012.

Liubartseva, S., Coppini, G., Lecci, R., and Clementi, E.: Tracking plastics in the Mediterranean: 2D Lagrangian model, Marine Pollution Bulletin, 129, 151–162, https://doi.org/10.1016/j.marpolbul.2018.02.019, 2018.

Lofty, J., Valero, D., Moreno-Rodenas, A., Belay, B. S., Wilson, C., Ouro, P., and Franca, M. J.: On the vertical structure of non-buoyant plastics in turbulent transport, Water Research, 254, 121306, https://doi.org/10.1016/j.watres.2024.121306, 2024.

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

Martin, J., Henrichs, T., Clark, B., Stanners, D., McInnes, G., McGlade, J., Petersen, J., Huntington, J., Vos, H., McAleavey, P., Uhel, R., Ribeiro, T., Gheorghe, A., Reichel, A.-D., Barbu, A., Jol, A., Barkman, A., Meiner, A., Lükewille, A., Mourelatou, A., Werner, B., Georgi, B., Kurnik, B., Romao, C., Adem, C., Gee, D., Clubb, D., Dejean, F., Vollmer, G., Aristei, G., Füssel, H.-M., Pereira, I., Weber, J.-L., Mortensen, L., Winograd, M., Erhard, M., Adams, M., Andersen, M. S., Asquith, M., Chrenko, Martin, J., Henrichs, T., Clark, B., Stanners, D., McInnes, G., McGlade, J., Petersen, J.-E., Huntington, J., Vos, H., McAleavey, P., Uhel, R., Ribeiro, T., Gheorghe, A., Reichel, A., Barbu, A.-D., Jol, A., Barkman, A., Meiner, A., Lükewille, A., Mourelatou, A., Werner, B., Georgi, B., Kurnik, B., Romao, C., Adem, C., Gee, D., Clubb, D. O., Dejean, F., Vollmer, G., Aristei, G., Füssel, H.-M., Pereira, I., Weber, J.-L., Mortensen, L., Winograd, M., Erhard, M., Adams, M., Skou Andersen, M., Asquith, M., Chrenko, M., Bock, N., Kazmierczyk, P., Jensen, P., Kristensen, P., Spyropoulou, R., Fernandez, R., Collins, R., Pignatelli, R., Speck, S., Isoard, S., Christiansen, T., Foltescu, V., and Laporte, V.: The European environment-state and outlook 2010: Synthesis, Publications Office of the European Union, https://doi.org/10.2800/45773, 2010.

Martínez-Gomariz, E., Russo, B., Gómez, M., and Plumed, A.: An approach to the modelling of stability of waste containers during urban flooding, Journal of Flood Risk Management, 13, e12558, https://doi.org/10.1111/jfr3.12558, 2020.

Mellink, Y. A. M., van Emmerik, T. H. M., and Mani, T.: Wind- and rain-driven macroplastic mobilization and transport on land, Scientific Reports, 14, 3898, https://doi.org/10.1038/s41598-024-53971-8, 2024.

Merritt, D. M. and Wohl, E. E.: Processes governing hydrochory along rivers: Hydraulics, hydrology, and dispersal phenology, Ecological Applications, 12, 1071–1087, https://doi.org/10.1890/1051-0761(2002)012[1071:PGHARH]2.0.CO;2, 2002.

Molazadeh, M., Calabro, G., Liu, F., Tassin, B., Rovelli, L., Lorke, A., Dris, R., and Vollertsen, J.: The role of turbulence in the deposition of intrinsically buoyant

MPs, Science of The Total Environment, 911, 168540, https://doi.org/10.1016/j.scitotenv.2023.168540, 2024.

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

Morales-Hernández, M., Murillo, J., and García-Navarro, P.: Diffusion-dispersion numerical discretization for solute transport in 2D transient shallow-flows, Environmental Fluid Mechanics, 9, 1217–1234, 2019.

Morales-Hernández, M., Sharif, M. B., Kalyanapu, A., Ghafoor, S., Dullo, T., Gangrade, S., Kao, S.-C., Norman, M., and Evans, K.: TRITON: A Multi-GPU open source 2D hydrodynamic flood model, Environmental Modelling & Software, 141, 105034, https://doi.org/10.1016/j.envsoft.2021.105034, 2021.

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

Nordam, T., Kristiansen, R., Nepstad, R., van Sebille, E., and Booth, A. M.: A comparison of Eulerian and Lagrangian methods for vertical particle transport in the water column, Geosci. Model Dev., 16, 5339–5363, https://doi.org/10.5194/gmd-16-5339-2023, 2023.

OECD: Financial Management of Flood Risk, OECD Publishing, Paris, https://doi.org/10.1787/9789264257689-en, 2016.

Olcina, J., Sauri, D., Hernández, M., and Ribas, A.: Flood policy in Spain: a review for the period 1983–2013, Disaster Prevention and Management, 25, 41–58, 2016.

Peeters, F. and Hofmann, H.: Length-Scale Dependence of Horizontal Dispersion in the Surface Water of Lakes, Limnology and Oceanography, 60, 1917–1934, 2015.

Persi, E., Petaccia, G., and Sibilla, S.: Large wood transport modelling by a coupled Eulerian–Lagrangian approach, Natural Hazards, 91, 59–74, 2018a.

Persi, E., Petaccia, G., Sibilla, S., Brufau, P., and García-Navarro, P.: Calibration of a dynamic Eulerian-lagrangian model for the computation of wood cylinders transport in shallow water flow, Journal of Hydroinformatics, 21, 164–179, https://doi.org/10.2166/hydro.2018.085, 2018b.

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

Portillo De Arbeloa, N. and Marzadri, A.: Modeling the transport of microplastics along river networks, Science of The Total Environment, 911, 168227, https://doi.org/10.1016/j.scitotenv.2023.168227, 2024.

Ripple, W. J., Wolf, C., Newsome, T. M., Barnard, P., and Moomaw, W. R.: World scientists’ warning of a climate emergency, BioScience, 70, 8–100, 2020.

Rutherford, J.: River Mixing, Wiley, Chichester, https://doi.org/10.1002/aheh.19950230614, 1994.

Sampson, C. C., Smith, A. M., Bates, P. D., Neal, J. C., Alfieri, L., and Freer, J. E.: A high-resolution global flood hazard model, Water Resources Research, 51, 7358–7381, https://doi.org/10.1002/2015WR016954, 2015.

Schreyers, L. J., van Emmerik, T. H., Huthoff, F., Collas, F. P., Wegman, C., Vriend, P., Boon, A., de Winter, W., Oswald, S. B.,
Schoor, M. M., Wallerstein, N., van der Ploeg, M., and Uijlenhoet, R.: River plastic transport and storage budget, Water Research, 121786, https://doi.org/10.1016/j.watres.2024.121786, 2024.

Soares-Frazão, S. and Zech, Y.: Dam-break flow through an idealised city, Journal of Hydraulic Research, 46, 648–658, https://doi.org/10.3826/jhr.2008.3164, 2008.

Sánchez-Salas, J., Flores, J., Jurado, E., Sáenz-Mata, J., Orozco-Figueroa, P., and Muro Pérez, G.: Hidrocoria en semillas de Agave victoriae-reginae T. Moore, especie en peligro de extinción: Morfología y anatomía como facilitadores de la hidrodispersión y germinación, Gayana. Botánica, 74, 251–261, https://doi.org/10.4067/S0717-66432017000200251, 2017.

Thielen, J., Bartholmes, J., Ramos, M.-H., and de Roo, A.: The European Flood Alert System – Part 1: Concept and development, Hydrol. Earth Syst. Sci., 13, 125–140, https://doi.org/10.5194/hess-13-125-2009, 2009.

Thomas, H. and Nisbet, T.: Modelling the hydraulic impact of reintroducing large woody debris into watercourses, Journal of Flood Risk Management, 5, 164–174, https://doi.org/10.1111/j.1753-318X.2012.01137.x, 2012.

United Nations: Sendai Framework for Disaster Risk Reduction 2015–2030, United Nations, https://www.undrr.org/publication/sendai-framework-disaster-risk-reduction-2015-2030 (last access: 9 October 2025), 2015.

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

Valero, D., Bayón, A., and Franca, M. J.: Urban Flood Drifters (UFDs): Onset of movement, Science of The Total Environment, 927, 171568, https://doi.org/10.1016/j.scitotenv.2024.171568, 2024.

Vallés, P., Morales-Hernández, M., García-Navarro, P., Roeber, V., and Caviedes-Voullième, D.: Enhancing Flood Analysis with a Lagrangian Transport Modeling and SERGHEI, in: Advances in Hydroinformatics—SimHydro 2023, Volume 1, 387–405, Springer, https://doi.org/10.1007/978-981-97-4072-7, 2023.

Vallés, P., Fernández-Pato, J., Morales-Hernández, M., Echeverribar, I., and García-Navarro, P.: A 2D shallow water flow model with 1D internal boundary condition for subgrid-scale topography, Advances in Water Resources, 189, 104716, https://doi.org/10.1016/j.advwatres.2024.104716, 2024.

Vallés, P., Caviedes-Voullième, D., and Morales-Hernández, M.: SERGHEI-LPT-RK, Zenodo [code], https://doi.org/10.5281/zenodo.14870918, 2025a.

Vallés, P., Caviedes-Voullième, D., and Morales-Hernández, M.: SERGHEI v2.1.0, Zenodo [code], https://doi.org/10.5281/zenodo.14871005, 2025b.

van Emmerik, T., Mellink, Y., Hauk, R., Waldschlaeger, K., and Schreyers, L.: Rivers as Plastic Reservoirs, Frontiers in Water, 3, 786936, https://doi.org/10.3389/frwa.2021.786936, 2022.

Wallemacq, P., Guha-Sapir, D., McClean, D., CRED, and UNISDR: The Human Cost of Weather-Related Disasters 1995–2015, Centre for Research on the Epidemiology of Disasters (CRED), Brussels, https://www.unisdr.org/files/46796_cop21weatherdisastersreport2015.pdf (last access: 9 October 2025), 2015.

Xia, J., Teo, F. Y., Lin, B., and Falconer, R. A.: Formula of incipient velocity for flooded vehicles, Natural Hazards, 58, 1–14, 2011.

Xia, X., Liang, Q., and Ming, X.: A full-scale fluvial flood modelling framework based on a high-performance integrated hydrodynamic modelling system (HiPIMS), Advances in Water Resources, https://doi.org/10.1016/j.advwatres.2019.103392, 2019.

Yang, H. and Foroutan, H.: Effects of near-bed turbulence on microplastics fate and transport in streams, Science of The Total Environment, 905, 167173, https://doi.org/10.1016/j.scitotenv.2023.167173, 2023.

Zamora, C. O. and Montagnini, F.: Seed Rain and Seed Dispersal Agents in Pure and Mixed Plantations of Native Trees and Abandoned Pastures at La Selva Biological Station, Costa Rica, Restoration Ecology, 15, 453–461, https://doi.org/10.1111/j.1526-100X.2007.00241.x, 2007.

Zhao, C., Fang, H., Ouro, P., Stoesser, T., and Dey, S.: Response of Bedload and Bedforms to Near-Bed Flow Structures, Journal of Hydraulic Engineering, 150, 04023060, https://doi.org/10.1061/JHEND8.HYENG-13618, 2024.