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

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 completeness of these tools, as object transport is a critical component of floods and can cause obstructions in hydraulic structures 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 debris accumulation near weirs, culverts and dams, highlighting their significant impact on flood risk (Thomas and Nisbet, 2012; De Cicco et al., 2015). Moreover, the transport of objects 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 others.

The transport of objects by fluid flows has been extensively studied in recent years (Braudrick and Grant, 2000; Xia et al., 2011; Martínez-Gomariz et al., 2020). These studies 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 transport 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 deposition, fragmentation and degradation, provided that these processes are explicitly implemented. On the other hand, the Eulerian method provides a different approach by treating debris as a concentration, offering greater computational efficiency 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 sometimes coastal scenarios with fixed dry–wet boundaries (Lebreton 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 feasible (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 contexts. 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 approach is to determine the temporal evolution of material by solving the Eulerian advection–dispersion equation, thereby obtaining the changes in debris concentration over time (Baharvand et al., 2023; Schreyers et al., 2024; Portillo De Arbeloa and Marzadri, 2024). While these models effectively capture 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 hydrodynamic frameworks (Jalón-Rojas et al., 2019; Pilechi et al., 2022), which provide more detailed flow information compared 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 Lagrangian 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 forecasting, real-time decision support, environmental impact assessments 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 water flow scenarios while incurring significantly lower computational 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 materials 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 transport (LPT) model, together with an accompanying analysis of the accuracy and computational efficiency of the numerical scheme and its implementation for particle transport. The accuracy and efficiency of an explicit Euler method are compared 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 intervals. Furthermore, particle dispersion is incorporated into the LPT model equations using a random-walk model because 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 Ecohydrology 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):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial h}{\partial t}}+\mathrm{\nabla }\left(h\mathit{v}\right)=r-i+e+S_w,\text{(2)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial hu}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left(u^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}}\right)+{\displaystyle \frac{\partial \left(huv\right)}{\partial y}}\\ & =-gh{S}_{\mathrm{f}x}-gh{\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial x}},\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial hv}{\partial t}}+{\displaystyle \frac{\partial }{\partial y}}\left(v^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}g{h}^{\mathrm{2}}\right)+{\displaystyle \frac{\partial \left(huv\right)}{\partial x}}\\ & =-gh{S}_{\mathrm{f}y}-gh{\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial y}},\end{array}\end{array}$$

where h is the water depth [L], $\mathit{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

$$\begin{array}{}\text{(4)}& S_{\mathrm{f}x}={\displaystyle \frac{n^{\mathrm{2}}u\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}}{h^{\mathrm{4}/\mathrm{3}}}},S_{\mathrm{f}y}={\displaystyle \frac{n^{\mathrm{2}}v\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}}{h^{\mathrm{4}/\mathrm{3}}}},\end{array}$$

where n is Manning's roughness coefficient [T L${}^{-\mathrm{1}/\mathrm{3}}$] (Arcement 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 ${\mathit{x}}_p=(x_p,y_p,z_p)$:

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{ll}\frac{\mathrm{d}{x}_p}{\mathrm{d}t}& =u\left({\mathit{x}}_p\right),\\ \frac{\mathrm{d}{y}_p}{\mathrm{d}t}& =v\left({\mathit{x}}_p\right),\\ z_p& =h\left({\mathit{x}}_p\right)+z_{\mathrm{b}}\left({\mathit{x}}_p\right).\end{array}\right.\end{array}$$

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:

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{ll}\frac{\mathrm{d}{x}_p}{\mathrm{d}t}& =u\left({\mathit{x}}_p\right)+u_{\text{disp}}\left({\mathit{x}}_p\right),\\ \frac{\mathrm{d}{y}_p}{\mathrm{d}t}& =v\left({\mathit{x}}_p\right)+v_{\text{disp}}\left({\mathit{x}}_p\right),\\ z_p& =h\left({\mathit{x}}_p\right)+z_{\mathrm{b}}\left({\mathit{x}}_p\right),\end{array}\right.\end{array}$$

where ${\mathit{v}}_{\text{disp}}=(u_{\text{disp}},v_{\text{disp}},\mathrm{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² 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. Therefore, the standard anisotropic diffusion model (Morales-Hernández et al., 2019) is utilized to determine the diffusivity variables, with the following expressions:

$$\begin{array}{}\text{(7)}& K_{hx}={\mathit{ϵ}}_L\left|u^{*}\right|h(x_p,y_p),K_{hy}={\mathit{ϵ}}_T\left|u^{*}\right|h(x_p,y_p),\end{array}$$

where ϵ_L, and ϵ_T are the longitudinal and transversal dispersion coefficients [L² T⁻¹] (Rutherford, 1994), respectively; and $\left|u^{*}\right|$ is the friction velocity:

$$\begin{array}{}\text{(8)}& \left|u^{*}\right|=n\sqrt{g{\displaystyle \frac{v_x^{\mathrm{2}}+v_y^{\mathrm{2}}}{\left(h\right(x_p,y_p){)}^{\mathrm{1}/\mathrm{3}}}}}.\end{array}$$

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ème 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:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{ll}{x}_p^{n+\mathrm{1}}& =x_p^n+\mathrm{\Delta }{t}_p^n{u}_i^n+\mathrm{\Delta }{t}_p^n{u}_{\text{disp},p}^n\\ & =x_p^n+\mathrm{\Delta }{t}_p^n{u}_i^n+R_{p,x}^n{\left(\mathrm{2}{\mathit{\sigma }}^{-\mathrm{1}}{K}_{hx,i}^n\mathrm{\Delta }{t}_p^n\right)}^{\mathrm{1}/\mathrm{2}},\\ y_p^{n+\mathrm{1}}& =y_p^n+\mathrm{\Delta }{t}_p^n{v}_i^n+\mathrm{\Delta }{t}_p^n{v}_{\text{disp},p}^n\\ & =y_p^n+\mathrm{\Delta }{t}_p^n{v}_i^n+R_{p,y}^n{\left(\mathrm{2}{\mathit{\sigma }}^{-\mathrm{1}}{K}_{hy,i}^n\mathrm{\Delta }{t}_p^n\right)}^{\mathrm{1}/\mathrm{2}},\\ z_p^{n+\mathrm{1}}& =h_i^n+z_{\mathrm{b},i},\end{array}\right.\end{array}$$

with the dispersion velocities defined by

$$\begin{array}{}\text{(10)}& u_{\text{disp},p}^n=\sqrt{{\displaystyle \frac{\mathrm{2}{K}_{hx,i}^n}{\mathit{\sigma }\mathrm{\Delta }{t}_p^n}}},v_{\text{disp},p}^n=\sqrt{{\displaystyle \frac{\mathrm{2}{K}_{hy,i}^n}{\mathit{\sigma }\mathrm{\Delta }{t}_p^n}}},\end{array}$$

where the Eulerian velocity field ${\mathit{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 ${\mathit{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 $\mathrm{\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 σ=1; $K_{hx,i}^n$ and $K_{hy,i}^n$ are the horizontal diffusivity in $\widehat{x}$ and $\widehat{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{array}{}\text{(11)}& \begin{array}{rl}{\mathit{x}}_p^{n+\mathrm{1}}& ={\mathit{x}}_p^n+{\displaystyle \frac{\mathrm{\Delta }{t}_p^n}{\mathrm{6}}}(\mathit{v}\left({\mathit{x}}_p,t^n\right)+\mathrm{2}\mathit{v}\left({\mathit{x}}_p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}a,t^{n+\mathrm{1}/\mathrm{2}}\right)\\ & +\mathrm{2}\mathit{v}\left({\mathit{x}}_p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}b,t^{n+\mathrm{1}/\mathrm{2}}\right)+\mathit{v}\left({\mathit{x}}_p+c,t^{n+\mathrm{1}}\right)),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& a=\mathrm{\Delta }{t}_p^n\mathit{v}\left({\mathit{x}}_p,t^n\right),b=\mathrm{\Delta }{t}_p^n\mathit{v}\left({\mathit{x}}_p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}a,t^{n+\mathrm{1}/\mathrm{2}}\right),\\ & c=\mathrm{\Delta }{t}_p^n\mathit{v}\left({\mathit{x}}_p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}b,t^{n+\mathrm{1}/\mathrm{2}}\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\mathit{v}}^{n+\mathrm{1}/\mathrm{2}}\left({\mathit{x}}_p\right)& ={\displaystyle \frac{\mathrm{5}}{\mathrm{16}}}{\mathit{v}}^{n+\mathrm{1}}\left({\mathit{x}}_p\right)+{\displaystyle \frac{\mathrm{15}}{\mathrm{16}}}{\mathit{v}}^n\left({\mathit{x}}_p\right)\\ & -{\displaystyle \frac{\mathrm{5}}{\mathrm{16}}}{\mathit{v}}^{n-\mathrm{1}}\left({\mathit{x}}_p\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{16}}}{\mathit{v}}^{n-\mathrm{2}}\left({\mathit{x}}_p\right).\end{array}\end{array}$$

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.

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).

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 x_p in the same position of the point $\mathit{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 $\mathit{q}=(q_x,q_y)$ are the lower-left and lower-right corners (respectively) of the cell, b is the distance between q and the intersection edge point and $\widehat{\mathit{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 $\widehat{\mathit{v}}$.

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.

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 Δt_p of particle p is given by

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }{t}_p\le {\displaystyle \frac{d}{\left|\left|{\mathit{v}}_i\right|\right|}}.\end{array}$$

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

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }{t}_p\le \mathrm{\Delta }t,\end{array}$$

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

$$\begin{array}{}\text{(16)}& \mathrm{\Delta }t=\sum _{m=\mathrm{1}}^M\mathrm{\Delta }{t}_{p,m},\end{array}$$

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 Δt. Since particles can travel through at most three cells, the number of subdivisions M satisfies $\mathrm{1}\le M\le \mathrm{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.

Figure 3 Three possible situations for a particle: traveling through 1 cell (a), 2 cells (b), or 3 cells (c).

The complete formulation for the particle time-step restrictions is

$$\begin{array}{}\text{(17)}& \mathrm{\Delta }{t}_{p,m}\le \text{min}\left({\displaystyle \frac{d}{\left|\left|{\mathbf{v}}_i\right|\right|}},\mathrm{\Delta }t-\sum _{m=\mathrm{1}}^{M-\mathrm{1}}\mathrm{\Delta }{t}_{p,m}\right),\end{array}$$

where i refers to the current cell containing the particle, which may change during the hydrodynamic time step Δ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 $\widehat{\mathit{v}}=\mathit{v}{\left|\left|\mathit{v}\right|\right|}^{-\mathrm{1}}$ 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):

   $$\begin{array}{}\text{(18)}& \mathit{p}+d\widehat{\mathit{v}}=\mathit{q}+b\widehat{e}.\end{array}$$

   Considering only the x-coordinate:

   $$\begin{array}{}\text{(19)}& p_x+d\widehat{v_x}=q_x,\end{array}$$

   and therefore, if $\widehat{v_x}$ is not null,

   $$\begin{array}{}\text{(20)}& d={\displaystyle \frac{q_x-p_x}{\widehat{v_x}}}={\displaystyle \frac{A_x+\mathrm{\Delta }x-p_x}{\widehat{v_x}}},\end{array}$$

   where Δ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 $\widehat{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):

   $$\begin{array}{}\text{(21)}& {\displaystyle }{d}_{\mathrm{2}}={\displaystyle \frac{A_x-p_x}{\widehat{v_x}}}\text{if and only if}\widehat{v_x}\ne \mathrm{0},\text{(22)}& {\displaystyle }{d}_{\mathrm{3}}={\displaystyle \frac{A_y+\mathrm{\Delta }x-p_y}{\widehat{v_y}}}\text{if and only if}\widehat{v_y}\ne \mathrm{0},\text{(23)}& {\displaystyle }{d}_{\mathrm{4}}={\displaystyle \frac{A_y-p_y}{\widehat{v_y}}}\text{if and only if}\widehat{v_y}\ne \mathrm{0}.\end{array}$$

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 nIter hydrodynamic time steps (from n to n+nIter), 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 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-Voulliè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.

Figure 4 Conceptual diagram for the SERGHEI modules and their implementation status.

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 hydrodynamic 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.

Figure 5 Schematic of the numerical calculation using the LPT module.

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.

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.

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

$$\begin{array}{}\text{(24)}& L_{\mathrm{1}}\left(t^n\right)=|{\mathit{x}}_p^n-\widehat{\mathit{x}}|,\end{array}$$

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

$$\begin{array}{}\text{(25)}& \text{MAE}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=\mathrm{0}}^N{L}_{\mathrm{1}}\left(t^n\right),\end{array}$$

where N is the number of time measurements.

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

$$\begin{array}{}\text{(26)}& L_{\mathrm{2}}\left(t^n\right)=({\mathit{x}}_p^n-\widehat{\mathit{x}}{)}^{\mathrm{2}},\end{array}$$

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

$$\begin{array}{}\text{(27)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{n=\mathrm{0}}^N{L}_{\mathrm{2}}\left(t^n\right)}.\end{array}$$

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 set-up 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 buildings. Finally, the differences between the methods' results are shown in Fig. 9, where the L₂-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.

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.

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⁵ particles in the dam break case.

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).

Figure 9 Temporal evolution of the L₂-norm per particle obtained from the difference between methods results.

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.

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).

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).

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.

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.

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.

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 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×10⁵ computational cells; 10⁴ 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 precipitation 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.

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.

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.

Figure 16 Temporal evolution of the total travel distance for both 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 insufficient 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 facilitate 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 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.