# A Robust Model for the Assessment of Oil Spill Hazards over Land and Water Bodies

Pablo Vallés, Sergio Martínez-Aranda, Reinaldo García, Pilar García-Navarro

► To cite this version:

Pablo Vallés, Sergio Martínez-Aranda, Reinaldo García, Pilar García-Navarro. A Robust Model for the Assessment of Oil Spill Hazards over Land and Water Bodies. Water, 2024, 16 (23), pp.3377. 10.3390/w16233377. hal-05317593

HAL Id: hal-05317593

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

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

Article

# A Robust Model for the Assessment of Oil Spill Hazards over Land and Water Bodies

Pablo Vallés¹,*, Sergio Martínez-Aranda ¹,*, Reinaldo García² and Pilar García-Navarro ¹

1 Fluid Dynamic Technologies TFD-I3A, Universidad de Zaragoza, 50018 Zaragoza, Spain

2 Hydronia LLC, Pembroke Pines, FL 33029, USA

* Correspondence: sermar@unizar.es

**Abstract:** Oil spills over land and water bodies are some of the most relevant hazards that should be considered when implementing oil production and transport projects. However, the development of robust, versatile, and efficient tools for carrying out this type of hazard assessment is a challenge for geophysical modellers due to the complexity of the oil flow over hybrid terrain-water surfaces. This work presents a versatile Eulerian approach to simulating the transport of an oil layer flowing over steep terrain that may also be dragged by an underlying water flow, i.e., rivers, lakes, oceans, etc., if it exists. The model allows for the seamless simulation of spills that start on land and eventually impact a water body in a single simulation step. The focus here is paid to the integration of the drag shear stresses between the layers, responsible for the oil spreading over a moving water surface. This drag term is solved using a non-iterative implicit method that allows for robust and efficient solutions even with high coupling between both layers. Two synthetic test cases are simulated to demonstrate the accuracy and robustness of the proposed model, obtaining results that validate the model's behaviour in high-coupling cases. Finally, the spreading hazard for a realistic oil production project is assessed. The results obtained verify the capability of the model to become a useful tool for oil spill forecasting over hybrid terrain-water surfaces.

**Keywords:** oil spill; surface shallow flows; drag shear stress; hybrid land-water surfaces; Eulerian approach

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

check for
updates

**Citation:** Vallés, P.; Martínez-Aranda, S.; García, R.; García-Navarro, P. A Robust Model for the Assessment of Oil Spill Hazards over Land and Water Bodies. *Water* **2024**, *16*, 3377. https://doi.org/10.3390/w16233377

Academic Editors: Xiang Yu, Yuke Wang and Yongqian Qu

Received: 14 October 2024

Revised: 18 November 2024

Accepted: 19 November 2024

Published: 24 November 2024

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

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Oil spills over land and water bodies are some of the most harmful hazards that should be considered when implementing development and industrial projects concerning oil transport. Although the number of oil spills has decreased globally decade after decade [1], these kinds of accidents continue to occur, and the transport of oil by oceans, rivers, and other water currents has increased in recent years, making the danger of a potential oil spill a real one [2]. In this sense, the necessity of robust, accurate, and efficient computational tools able to simulate potential hazards in multiple plausible scenarios is justified.

The vast majority of the simulation models for oil spills over water are based on a Lagrangian approach, focused on using particle tracking algorithms to represent the oil slick fate [3]. The MIKE 21/3 Oil Spill model [4] simulates oil spills on water bodies, including key weathering processes. However, it is not suitable for hybrid land-water spreading scenarios and lacks GPU acceleration for large-scale simulations. The Automated Cellular Lagrangian Network (ACLN) model for oil spills [5] allows the tracking of oil spreading in water environments. Cell Automata models use raster-type data structures and can be integrated with a GIS easily to provide an environmental database and display the results, but they use an extremely simplified formulation for the oil dynamics. The SaTCN model [6] combines a self-attention network with a temporal convolutional network to obtain a corrected estimation of historical wind fields and perform accurate backtracking of oil spills in a timely manner. The OILMAP model [7] is an oil spill response model
which provides algorithms for oil spreading, evaporation, and emulsification in water and shoreline domains. Although Lagrangian approaches to oil spills provide a high computational efficiency, they often overlook the interaction between the oil layer and the underlying layers, such as water bodies or the terrain. This drawback limits their ability to capture the intricate processes involved in the propagation of oil spills in hybrid land-aquatic areas [8].

For this reason, several Eulerian models [9] and Eulerian-Lagrangian models [10,11] have been developed to predict oil spills over water bodies. In this approach, the problem of an oil spill over water is usually modelled using two-layer mathematical formulations, which require complex numerical schemes and high computational effort [12,13]. A fully Eulerian numerical model for the simulation of such environmentally significant disasters was implemented by means of a two-dimensional (2D), two-layer shallow flow model in [14]. The authors considered a very thin oil layer over a thicker water layer in order to neglect the pressure term that the oil layer exerts on the water. Despite all of these advances, addressing robust and efficient high-resolution models that are able to deal with oil spills over hybrid land-water surfaces is still a challenge.

The present work is devoted to an even simpler Eulerian approach where the oil layer runs over land but can also be dragged by a underlying water flow if it exists (the bottom shear layer, referred to as BSL). The focus here is paid to the relative drag shear stresses between layers, responsible for the features of the oil over a moving volume of water. The numerical solution adopted is based on a finite volume (FV) upwind scheme with a Roe solver for the oil layer, solving its thickness, density, and horizontal velocity. Special care is taken to numerically integrate the drag terms at the oil-water interface, and a novel implicit procedure is proposed. The objective of the model is to be efficient, accurate, and robust in all situations, made possible by the previously calculated water flow movement, including the presence of wet-dry fronts in the oil layer moving over hybrid terrain-water surfaces.

The rest of this paper is structured as follows: the governing equations used for the description of the oil layer's movement are first presented in Section 2. Special emphasis is devoted to the formulation of the shear stress in the lower part of the oil column, including the friction force with the terrain if necessary and the drag shear stress induced by the water below. Then, in Section 3, the finite volume (FV) method deployed to solve the partial differential equations is outlined, with mention of the technique used to properly incorporate the drag shear stresses. Section 4 is devoted to the numerical results. The proposed method is first applied to two synthetic tests specifically designed to demonstrate the capabilities of the proposed approach. Then, the model is tested for the assessment of spreading hazards for a realistic case of an oil spill over land-water. Finally, the main conclusions are drawn in Section 5.

The model proposed here has been included as a new development within the Hydronia LLC software OilFlow2D (https://www.hydronia.com/oilflow2d, accessed on 1 September 2024), and some of the capabilities of the parent software were used in the testing cases. It is worth mentioning that oil processes, such as evaporation, dispersion, or emulsification, are not included in this work. This work is focused on the particular problem of the numerical discretisation of the drag friction forces, ensuring the stability and performance of the numerical model.

## 2. System of Equations

In the depth-averaged equations for shallow flows on environmental surfaces, the pressure is assumed to be hydrostatic, and the contributions of lateral stresses $T_{ij}$ and dispersive terms $D_{ij}$ are considered negligible. The resulting system of 2D depth-averaged equations, focusing on the governing equations of the motion of a thin oil layer on water or land, can be expressed using the following notation [15]:

$$ \frac{\partial(\bar{\rho}h)}{\partial t} + \frac{\partial(\bar{\rho}h\bar{u})}{\partial x} + \frac{\partial(\rho h \bar{v})}{\partial y} = 0 \quad (1) $$
$$ \frac{\partial(\bar{\rho}h\bar{u})}{\partial t} + \frac{\partial}{\partial x} \left( \bar{\rho}h\bar{u}^2 + \frac{1}{2}g\bar{\rho}h^2 \right) + \frac{\partial(\bar{\rho}h\bar{u}\bar{v})}{\partial y} = -\bar{\rho}g\bar{h} \frac{\partial z_b}{\partial x} - \tau_x + \tau_{slx} \quad (2) $$

$$ \frac{\partial(\bar{\rho}h\bar{v})}{\partial t} + \frac{\partial(\bar{\rho}h\bar{u}\bar{v})}{\partial x} + \frac{\partial}{\partial y} \left( \bar{\rho}h\bar{v}^2 + \frac{1}{2}g\bar{\rho}h^2 \right) = -\bar{\rho}g\bar{h} \frac{\partial z_b}{\partial y} - \tau_y + \tau_{sly} \quad (3) $$

where $\bar{\rho}$ is the depth-averaged oil density; $h$ denotes the oil depth; $\bar{\rho}h\bar{u}$ and $\bar{\rho}h\bar{v}$ are the unit mass flow rates along the x- and y-coordinates, respectively, where $\bar{u} = (\bar{u}, \bar{v})$ is the depth-averaged oil flow velocity; $z_b$ is the lower interface elevation, i.e., the terrain or the free surface of the bottom water layer; and $g$ is the gravitational acceleration.

The terms $\tau_x$ and $\tau_y$ are the horizontal components of the friction stress between the oil flow and the terrain surface if the oil spill is over land. This frictional stress can be expressed as $\tau = (\tau_x, \tau_y) = \tau \mathbf{n}_u$, where $\tau$ is the frictional loss and $\mathbf{n}_u = (n_{ux}, n_{uy}) = \bar{u} / |\bar{u}|$ is the unit relative velocity direction vector of the oil layer, where $|\bar{u}| = \sqrt{\bar{u}^2 + \bar{v}^2}$ is the oil velocity modulus. The basal frictional loss in oil-over-land depth-averaged models applied to Earth surface processes is assumed mainly viscous and governed by the simplified Bingham formula [16] as

$$ \tau = \frac{3}{2} \tau_y + 3\mu \frac{|\bar{u}|}{h} \quad (4) $$

where $\tau_y$ is the yield strength and $\mu$ is the oil viscosity.

The terms $\tau_{slx}$, $\tau_{sly}$ are the components of the drag shear stress at the lower surface due to the presence of an underlying water layer with a 2D velocity field, i.e., rivers, lakes, oceans, and other water bodies. This drag shear stress can be expressed as $\tau_{sl} = (\tau_{slx}, \tau_{sly}) = \tau_{sl} \mathbf{n}_r$, where $\tau_{sl}$ is the drag stress modulus and $\mathbf{n}_r = (n_{rx}, n_{ry}) = \bar{u} / |\bar{u}|$ is the direction of the relative velocity between the oil layer and the bottom shear layer $\mathbf{u}_r = \mathbf{u}_{sl} - \bar{u}$. Depending on the assumptions, the drag stress modulus must be expressed using different formulations. In this work, a turbulent formulation is used to define the drag stress modulus $\tau_{sl}$ as

$$ \tau_{sl} = \rho_{sl} C_{sl} |\mathbf{u}_r|^2 \quad (5) $$

where $\rho_{sl}$ is the density of the drag layer and $C_{sl}$ is the friction coefficient between the oil layer and the drag layer.

Equations (1)–(3) can be expressed in vector form as

$$ \frac{\partial \mathbf{U}}{\partial t} + \nabla \cdot \mathbf{E}(\mathbf{U}) = \mathbf{S}_b(\mathbf{U}) + \mathbf{S}_{sl}(\mathbf{U}) \quad (6) $$

where $\mathbf{U}$ is the vector of conservative variables:

$$ \mathbf{U} = (\bar{\rho}h, \bar{\rho}h\bar{u}, \bar{\rho}h\bar{v})^T \quad (7) $$

and $\mathbf{E}(\mathbf{U}) = (\mathbf{F}(\mathbf{U}), \mathbf{G}(\mathbf{U}))$ are the convective fluxes in the x- and y-directions, respectively:

$$ \mathbf{F}(\mathbf{U}) = \begin{pmatrix} \bar{\rho}h\bar{u} \\ \bar{\rho}h\bar{u}^2 + \frac{1}{2}g\bar{\rho}h^2 \\ \bar{\rho}h\bar{u}\bar{v} \end{pmatrix} \qquad \mathbf{G}(\mathbf{U}) = \begin{pmatrix} \bar{\rho}h\bar{v} \\ \bar{\rho}h\bar{u}\bar{v} \\ \bar{\rho}h\bar{v}^2 + \frac{1}{2}g\bar{\rho}h^2 \end{pmatrix} \quad (8) $$

The vector $\mathbf{S}_b(\mathbf{U})$ is defined by the momentum source terms associated with the lower surface slope and to the terrain friction loss, if it exists, whereas $\mathbf{S}_{sl}(\mathbf{U})$ is defined by the momentum source terms associated with the drag shear stress:

$$ \mathbf{S}_b(\mathbf{U}) = \begin{pmatrix} 0 \\ -\bar{\rho}g\bar{h} \frac{\partial z_b}{\partial x} - \tau n_{ux} \\ -\bar{\rho}g\bar{h} \frac{\partial z_b}{\partial y} - \tau n_{uy} \end{pmatrix} \qquad \mathbf{S}_{sl}(\mathbf{U}) = \begin{pmatrix} 0 \\ \tau_{sl} n_{rx} \\ \tau_{sl} n_{ry} \end{pmatrix} \quad (9) $$
### 3. Numerical Methods

To solve the system of Equation (6), a numerical method is necessary because it is a system without an analytical solution. In this work, an explicit upwind finite volume numerical scheme based on Roe’s solver is used to solve it [17,18]. The spatial domain is divided into triangular computational cells using a fixed mesh in time and the system (6) is integrated in each mesh cell $\Omega_i$. It is important to remark that, from this section onwards, averaging of the conservative variables $\bar{\rho}$, $\bar{u}$, and $\bar{v}$ is assumed and, therefore, overbar notation is no longer used for these variables. For more details on the numerical scheme, some references are included [15,17,19].

The expression to update the conservative variables from the current time level $t^n$ to the next time level $t^{n+1}$ is given by

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{A_i} \sum_{k=1}^{\text{NE}} \mathbf{R}_k^{-1} \mathbf{F}_k^{\downarrow} l_k + \frac{\Delta t}{A_i} \int_{\Omega_i} \mathbf{S}_{\text{sl}}(\mathbf{U}) \, d\Omega \quad (10) $$

where $\Delta t = t^{n+1} - t^n$ is the time step, $A_i$ is the discrete cell area, NE is the number of edges for the $i$ cell, $\mathbf{R}_k$ is the rotation matrix, $l_k$ is the length of the $k$th edge, and $\mathbf{F}_k^{\downarrow}$ is the explicit intercell flux vector at the $k$th cell edge, including the conservative fluxes, the bed-pressure, and the terrain friction momentum contributions [17].

The intercell flux vector $\mathbf{F}_k^{\downarrow} = (F_1, F_2, F_3)^{\downarrow}$ at the $k$th cell edge, separating the left $i$ cell and the right $j$ cell, is computed here as

$$ \mathbf{F}_k^{\downarrow} = \mathbf{F}(\mathbf{R}_k \mathbf{U}_i^n) + \sum_{m-} [(\tilde{\lambda}_m \tilde{\alpha}_m - \tilde{\beta}_m) \tilde{e}_m]_k^n \quad (11) $$

where $\tilde{\lambda}_{m,k}$ and $(\tilde{e}_m)_k$ are the eigenvalues and eigenvectors of the flux Jacobian matrix, respectively; $\tilde{\alpha}_{m,k}$ denotes the wave strengths; and $\tilde{\beta}_{m,k}$ are the source strengths which include the integrated bed-pressure and terrain friction contributions through the cell edge [17]. The subscript $m-$ under the sums indicate that only the waves travelling inward of the $i$ cell are considered, leading to the upwind computation of the flux at the edge.

It is important to note that the time integration is explicit, so the allowed value of the time step size $\Delta t$ must be limited to ensure stability. Here, the time step is dynamically controlled using the characteristic eigenvalues of the convective flux matrix and imposing the CFL $\le 1$ condition for triangular meshes, as described in [17]. This is an important drawback of explicit schemes. The present work is concerned with the proper discretisation of the drag shear stress term that enables stable computation, avoiding further time step reduction.

As can be seen in expression (10), all terms are discrete and can be obtained from the numerical scheme except the term accounting for the drag shear stress, which must be solved independently from the rest of the source terms. In this work, an implicit centred method for the numerical resolution of the cell-centred integral $\int_{\Omega_i} \mathbf{S}_{\text{sl}}(\mathbf{U}) \, d\Omega$ in (10) is presented.

#### Implicit Centred Method

Using an implicit centred method, the drag layer term in (10) can be expressed at time $t = t^{n+1}$ for each $i$ cell as

$$ \int_{\Omega_i} \mathbf{S}_{\text{sl}}(\mathbf{U}(x, y, t^{n+1})) \, d\Omega = A_i \mathbf{S}_{\text{sl}}(\mathbf{U}_i^{n+1}) \quad (12) $$

and the updated expression (10) can be written as

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{A_i} \sum_{k=1}^{\text{NE}} \mathbf{R}_k^{-1} \mathbf{F}_k^{\downarrow} l_k + \Delta t \mathbf{S}_{\text{sl}}(\mathbf{U}_i^{n+1}) \quad (13) $$
and in the flux term, the bed slope and terrain friction source terms, as well as the drag shear stress term $S_{sl}(U_i^{n+1})$, are computed separately as

$$S_{sl}(U_i^{n+1}) = \begin{pmatrix} 0 \\ \tau_{sl} n_{rx} \\ \tau_{sl} n_{ry} \end{pmatrix}_i^{n+1} \qquad (14)$$

From (13), the expressions for the provisional updating of the conservative mass rates $(\rho hu)_i^*$ and $(\rho hv)_i^*$ along the x- and y-coordinates, respectively, can be obtained as

$$(\rho hu)_i^* = (\rho hu)_i^n - \frac{\Delta t}{A_i} \sum_{k=1}^{NE} (F_2 n_x - F_3 n_y)_k^{\downarrow-} l_k \qquad (15)$$

$$(\rho hv)_i^* = (\rho hv)_i^n - \frac{\Delta t}{A_i} \sum_{k=1}^{NE} (F_2 n_y + F_3 n_x)_k^{\downarrow-} l_k$$

Using these terms of provisional mass rates (15) and operating the momentum equations in (13), it is possible to obtain a final updating expression of all the conservative variables, which can be written as

$$(\rho h)_i^{n+1} = (\rho h)_i^n - \frac{\Delta t}{A_i} \sum_{k=1}^{NE} (F_1)_k^{\downarrow-} l_k \qquad (16a)$$

$$(\rho hu)_i^{n+1} = \frac{(\rho hu)_i^* + B_i^* (\rho h)_i^{n+1} (u_{sl})_i^{n+1}}{1 + B_i^*} \qquad (16b)$$

$$(\rho hv)_i^{n+1} = \frac{(\rho hv)_i^* + B_i^* (\rho h)_i^{n+1} (v_{sl})_i^{n+1}}{1 + B_i^*} \qquad (16c)$$

where the factor $B_i^{n+1}$ is obtained using (5):

$$B_i^* = \frac{\rho_{sl} C_{sl} |\mathbf{u}_r|_i^n \Delta t}{(\rho h)_i^{n+1}} > 0 \qquad (17)$$

where $|\mathbf{u}_r|_i^n$ is the relative velocity modulus obtained using the value of the local variables at time $t^n$ and $(\rho h)_i^{n+1}$ is the oil mass in the cell at time $t^{n+1}$. It is worth mentioning that this proposed non-iterative implicit integration avoids the necessity of imposing additional time step restrictions to ensure the stability of the numerical solution, even for large drag shear stresses. Therefore, the condition CFL ≈ 1 can be imposed by default.

# 4. Results

This section presents the numerical results obtained for two synthetic tests specifically designed for the work, as well as for a realistic oil spill case, their interpretation, and the main conclusions that can be drawn.

## 4.1. Synthetic Oil Spill in a Straight Channel with Parabolic Velocity

This is a synthetic test for checking the behaviour of the proposed model in controlled conditions. A wide, flat, and straight rectangular channel of length $L = 1000$ m and width $W = 100$ m is considered. A steady velocity field $u_{sl} = (u_{sl}(y), 0)$ following a parabolic function $u_{sl}(y) = u_{sl}^{max} - (2y/W)^2$ along the channel cross-section is assumed to act as a bottom shear layer (BSL) for the oil spill. Therefore, the BSL velocity is null at the channel sides $u_{sl}(y = -50)$ m) = $u_{sl}(y = 50)$ m) = 0 and maximum at the channel centre $u_{sl}(y = 0)$ m) = $u_{sl}^{max} = 1$ m/s. The BSL velocity field is shown in Figure 1.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Bottom shear layer velocity field $\mathbf{u}_{sl} = (u_{sl}(y), 0)$.</figcaption></figure>

At the inlet section $x = 0$ m, a constant depth $h_{in} = 0.5$ m is imposed for the oil inflow with density $\rho = 860$ kg/m$^3$, viscosity $\mu = 0.001$ Pa·s, and yield strength $\tau_y = 10$ Pa. The proposed BSL model is tested for a bottom shear layer density $\rho_{sl} = 1000$ kg/m$^3$ and for drag coefficients varying from $C_{sl} = 10^{-4}$ to $C_{sl} = 10^2$. The computational domain is discretised by almost $2 \cdot 10^4$ triangular cells with a $10$ m$^2$ average area. The time step is limited using CFL = 0.95 and the simulation evolves until a steady solution for the oil flow is reached.

Figure 2 shows the oil depth along the longitudinal axis at the centre of the channel for all the drag coefficients $C_{sl}$ tested. The dashed black line shows the solution for the case with the BSL drag disabled. It is worth noting that, for this case, the frictional resistance $\tau_b$ is lower than the pressure force at the inlet section and does not allow the flow to develop. Hence, the oil layer stops at $x \approx 110$ m if the BSL is removed. When the bottom shear drag is enabled, the model is able to compute stable steady solutions with almost uniform oil depths along the channels, regardless of the order of the drag coefficient $C_{sl}$ considered. The uniform oil layer decreases as the drag coefficient increases, ranging from $h \approx 0.48$ m for $C_{sl} = 10^{-4}$ to $h \approx 0.46$ m for $C_{sl} = 10^2$.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Depth of the steady flow at the longitudinal profile $y = 0$ m with different values of the drag coefficient $C_{sl}$.</figcaption></figure>

Moreover, Figure 3 depicts the temporal evolution of the channel outflow for all the drag coefficients $C_{sl}$ tested. Note that the model is always able to compute a steady solution for the outlet discharge, regardless of the order of the drag coefficient $C_{sl}$ considered. All the simulations converge to approximately the same final discharge $Q_{out} \approx 31$ m$^3$/s without numerical oscillations. This fact shows the robustness of the proposed model, allowing highly coupled relative displacements to be computed between adjacent layers, maintaining the numerical stability. However, the convergence is slower as the drag coefficient is $C_{sl}$ decreases.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Outlet discharge evolution with different values of the drag coefficient $C_{sl}$.</figcaption></figure>

Contrarily, the coupling between the BSL velocity field and the oil layer velocity becomes stronger as the drag coefficient $C_{sl}$ is increased. Figure 4 shows the velocity modulus $|\mathbf{u}|$ of the oil layer once the steady state is reached for three of the drag coefficients $C_{sl}$ tested. It is worth noting that the higher the drag coefficient, the faster the adaptation of the oil velocity field $|\mathbf{u}|$ to the drag shear layer velocity $|\mathbf{u}_{sl}|$. For $C_{sl} \ge 1$, the coupling of the oil layer and the BSL is almost immediate.

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Velocity modulus field of the steady flow with different values of the drag coefficient $C_{sl}$.</figcaption></figure>

The enhanced layer coupling with the higher drag coefficients $C_{sl}$ is clearly shown by the cross-sectional velocity profile at the end of the channel. Figure 5 shows the depth-averaged longitudinal velocity $u$ of the oil layer along the $y$-coordinate at $x = 900$ m for all the drag coefficients $C_{sl}$ tested. Also, the BSL longitudinal velocity $u_{sl}$ profile is plotted with a dashed black line for comparison. For the lower drag coefficients $C_{sl}$, a marked slip appears between the BSL and the oil layer, leading to an almost uniform velocity profile in the oil layer at the end of the channel. Nevertheless, the higher the drag coefficient $C_{sl}$, the stronger the coupling between the BSL and the oil layer. For drag coefficients larger than $C_{sl} = 1$, the velocity profile of the oil layer at the end of the channel shows almost full agreement with the BSL velocity.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Velocity modulus of the steady flow at the cross-section $x = 900$ m with different values of the drag coefficient $C_{sl}$.</figcaption></figure>

## 4.2. Synthetic Oil Spill on a Rotating Bay

This synthetic test is designed to demonstrate the capability of the proposed model to address an oil spill that transits from land to water. The test geometry is a two-dimensional bay 1000 m long along the $x-$ and $y-$ coordinates, and a concentric terrain slope with elevations ranging from $z_b \in [0, 10]$ m as in Figure 6a. We consider a water bottom shear layer (BSL) with a constant and uniform free surface elevation at level 7.5 m. The BSL has a uniform velocity modulus $|u_{sl}| = 0.1$ m/s and the velocity direction rotates during 24 h, varying from the west to west directions progressively (Figure 6b). The oil inflow occurs at the left side of the domain, with a constant inlet discharge $Q_{in} = 0.2$ m$^3$/s.

The domain is discretised using an unstructured triangular mesh with 29,614 cells of a 50 m$^2$ median area. The simulated temporal window is 24 h and the condition CFL = 0.95 is set for the dynamical time step limitation. The west side of the domain is an inlet boundary for the oil layer inflow, whereas the north, east, and south sides are left as open outlet boundaries. The oil inflow at the west boundary is considered with density $\rho = 860$ kg/m$^3$ and viscosity $\mu = 0.1$ Pa·s. The water BSL is $\rho_{sl} = 1000$ kg/m$^3$ and characterised by two different drag coefficients $C_{sl} = 10^{-4}$ and $C_{sl} = 10^{-2}$ that have been tested for analysing the coupling of the rotating BSL and the oil layer.

<figure><img src="image_10.png" /><figcaption>(a) Terrain elevation and BSL surface level.</figcaption></figure>

**Figure 6.** Cont.
<figure><img src="image_11.png" /><figcaption>(b) Evolution of the BSL velocity direction during the simulated period.
**Figure 6.** Terrain and bottom shear layer setup.</figcaption></figure>

Figure 7 shows the oil layer depth (colour map) and velocity (vectors) at different times for a drag coefficient (left column) $C_{sl} = 10^{-4}$ and (right column) $C_{sl} = 10^{-2}$. For $x \in [0, 200]$, the oil moves over the original terrain following the higher bed slope direction and both simulations show a similar solution for all times. However, for $x \ge 200$ m, the oil layer reaches the rotating BSL and the oil velocity is influenced by the direction of the BSL velocity. The model allows the oil layer velocity direction to change following the rotating BSL, regardless of the drag coefficient $C_{sl}$ set.

<figure><img src="image_12.png" /><figcaption>**Figure 7.** *Cont.*</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Oil layer depth and velocity evolution for (left column) $C_{sl} = 10^{-4}$ and (right column) $C_{sl} = 10^{-2}$. From top to bottom: $t = 3$ h, $t = 9$ h, and $t = 15$ h.</figcaption></figure>

Furthermore, as the drag coefficient increases from $C_{sl} = 10^{-4}$ to $C_{sl} = 10^{-2}$, the coupling between the BSL and the oil layer is enhanced and the spreading of the oil within the bay is limited. Note that the oil layer is fully enclosed within the bay for the last times due to the BSL velocity direction for both drag coefficients $C_{sl}$ tested. This result demonstrates the capability of the model to include external boundary shear stresses into the 2D conservative oil flow system by means of a robust and accurate procedure.

The proposed implicit algorithm is also able to deal with wet-dry fronts for the oil layer without requiring additional time step limitations over those imposed by the CFL condition on the oil layer flow. Figure 8 shows the time step evolution during both simulations. The case with BSL disabled is also plotted for comparison. Note that the proposed procedure to include boundary shear stresses allows a similar time step limitation to be maintained without affecting the model stability, regardless of the drag coefficient $C_{sl}$ selected.
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Time step during the simulation.</figcaption></figure>

### 4.3. *Spill Hazard Assessment in Albert Lake (Uganda)*

This test aims to illustrate the capability of the proposed model to become a useful tool for the assessment of oil spill hazards, environmental impacts, and contingency plans, associated with realistic projects in hybrid land-water environments. In this sense, this section does not pretend to predict actual affected areas. As it is demonstrated below, the model provides a methodology that simplifies the performance of complex oil spills that start in land and eventually reach a water body. For the test, we consider the Tilenga Project, which is a component of a wider oil and gas infrastructure currently under development in the Albert Lake region (Uganda). All the data used in the following analysis result from reasonable estimations or have been obtained from web-information-based IA tools, and have not used proprietary information.

The Tilenga Project consists of 34 well pads, each containing up to 22 wells, connected by a network of buried pipelines that will transport the produced oil to a central processing facility. The expected oil production is estimated around 190 Mstb (3 × 10⁴ m³/day approximately) during the first six operation years [20]. The Tilenga Project area is situated close to the Victoria Nile estuary, one of Africa's most important regions for biodiversity and classified as a Critical Habitat by the International Finance Corporation (IFC). Figure 9 shows the Tilenga Project location, the well pad location (from Google satellite images), and the Victoria Nile-Albert Lake-White Nile hydrodynamic velocity with a steady inflow $Q_{Victoria-Nile}$ = 1700 m³/s at the upstream section of the Victoria Nile. The velocity modulus is lower than 1 m/s in most of the zones of the Victoria Nile and decreases considerably in the Albert lake and White Nile. This hydrodynamic steady state is considered the bottom shear layer 2D velocity field for the oil spill hazard assessment.

In this area, the project will implement three well fields located at the north (NF), at the southwest (SWF), and at the southeast (SEF) of the Victoria Nile estuary (see Figure 9). The characteristics of the oil in each well field are summarised in Table 1, including estimations of the oil type, the in situ viscosity, and the wax appearance temperature. The expected daily oil production for each well pad has been estimated by averaging the expected production of the whole Tilenga Project between the number of active well pads.

<table><caption>**Table 1.** Main characteristics of the oil and well pads for each field.</caption><thead><tr><th></th><th>North Field<br>NF</th><th>Southwest Field<br>SWF</th><th>Southeast<br>SEF</th></tr></thead><tbody><tr><td>Oil type</td><td>Heavy</td><td>Waxy</td><td>Waxy</td></tr><tr><td>Well pad production (m³/day)</td><td>~900</td><td>~900</td><td>~900</td></tr><tr><td>In situ viscosity (Pa·s)</td><td>~0.150</td><td>~0.005</td><td>~0.025</td></tr><tr><td>Wax appearance temperature (°C)</td><td>~30</td><td>~60</td><td>~45</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>**Figure 9.** Tilenga Project location and well pads, together with the hydrodynamic steady velocity field of the Victoria-White Nile rivers.</figcaption></figure>

For the spill hazard assessment, we select from each field the well pads closer to the Victoria Nile estuary: NF-01, NF-02, and NF-03 from the north field; SWF-01 and SWF-02 from the southwest field; and SEF-01 from the southeast field. Figure 10 depicts a zoom from the location of the selected well pads. Usually, the blowout hazard is assessed by considering spill scenarios with constant blowout rates over a 30-day duration. We use the generic spill scenarios and oil properties for estimating the spreading of the oil without considering mitigation measures. The BSL accounts for the hydrodynamic steady velocity computed previously for the Victoria Nile-Albert Lake-White Nile system and is modelled with water density $ρ_{sl} = 1000 \text{ kg/m}^3$ and a low drag coefficient $C_{sl} = 10^{-4}$. The spatial domain is discretised using an unstructured mesh with almost $1.5 \times 10^5$ triangular cells, with a 160 m maximum cell side and refined to a 50 m cell size in high-velocity zones. The time step is limited using CFL = 0.95 and the simulated temporal window is 30 days.

<figure><img src="image_16.png" /><figcaption>**Figure 10.** Selected well pads for the spill hazard assessment.</figcaption></figure>

Figure 11 depicts the oil spreading over a 30-day spill duration for the six well pads considered. Except for well pad NF-02, which is contained overland, the rest of the well pads show spills that reach the Victoria Nile. Once the oil layer makes contact with the
water current, the BSL model activates and the oil spreads downstream along the Albert Lake and the White Nile. In real applications, the results of the model similar to the ones presented herein could be used to support mitigation plans to contain the oil spreading and minimise environmental impacts.

<figure><img src="image_17.png" /><figcaption>**Figure 11.** Oil layer depth over a 30-day spill duration for blowout events in well pad (top-left) NF-01, (top-right) NF-02, (centre-left) NF-03, (centre-right) SWF-01, (bottom-left) SWF-02, and (bottom-right) SEF-01.</figcaption></figure>

The travelled distances along the Victoria Nile–Albert Lake–White Nile and the computational time required by the simulations are summarised in Table 2. Simulations were performed with OilFlow2D-GPU on a device with an NVIDIA GeForce GTX 1080-Ti (NVIDIA Corporation, Santa Clara, CA, USA, 2017). The maximum spill over a 30-day duration is obtained for the SWF-01 well pad, with a spreading distance about 32 km. The computational times are almost similar for all the spill scenarios and always lower than 20 min, even with this medium-grade GPU device, demonstrating the high efficiency of the proposed model. Finally, it is worth noting that the present application simplifies extremely complex physical–chemical processes and does not consider surface features such as vegetation and soil properties. Therefore, these results could overestimate the real impact of blowout events.
<table><caption>**Table 2.** Spreading distances and computational times for the blowout spills over a 30-day duration.</caption><thead><tr><th>Well Pads</th><th>Spreading Distances (km)</th><th>Comput. Time (min's'')</th></tr></thead><tbody><tr><td>NF-01</td><td>26.7</td><td>16' 00"</td></tr><tr><td>NF-02</td><td>2.2</td><td>17' 18"</td></tr><tr><td>NF-03</td><td>14.2</td><td>15' 09"</td></tr><tr><td>SWF-01</td><td>31.8</td><td>15' 25"</td></tr><tr><td>SWF-02</td><td>26.9</td><td>14' 00"</td></tr><tr><td>SEF-01</td><td>30.7</td><td>19' 12"</td></tr></tbody></table>

## 5. Conclusions

An Eulerian model has been used to model oil spills on land and water, forming a very thin layer and focusing on the drag forces exerted by the solid or liquid surfaces. A depth-averaged formulation of the shallow water type has been followed. The resulting hyperbolic system facilitates the development of robust numerical models capable of handling large-scale geophysical events where the water flow below the oil layer is known.

A series of cases have been performed to demonstrate the ability of the model to reproduce the correct oil transport. A first test case is simulated to verify the accuracy of the model under controlled conditions, computing highly coupled relative displacements between adjacent layers while maintaining numerical stability. The second test case is simulated to obtain a demonstration of the model's ability to simulate an oil spill transiting from land to water. The proposed model can deal with wet-dry fronts for the hydrocarbon layer without requiring additional time step constraints to those imposed by the CFL condition on the hydrocarbon layer flow as can be observed in the results. Finally, the spreading hazard for a realistic oil production project is assessed. The obtained results verify the capability of the model to become a useful tool for oil spill forecasting over hybrid terrain-water surfaces.

In summary, the proposed model provides promising results in the field of oil spill modelling, and could be an effective tool to prevent major damage in the event of a spill. In this sense, the model has been implemented as a new development for the Hydronia LLC software OilFlow2D (https://www.hydronia.com/oilflow2d, accessed on 1 September 2024). The high computational efficiency achieved thanks to the GPU acceleration and the lack of additional time step limitations due to the implicit integration of the drag shear stresses allows an affordable computational effort even in large-scale long-term scenarios. However, the possibility of developing a numerical resolution for the drag layer terms by other numerical integration methods remains as future work, to improve accuracy, robustness, and efficiency of the computational tools available for these geophysical processes.

**Author Contributions:** Conceptualization, R.G. and P.V.; methodology, P.V. and S.M.-A.; software, S.M.-A.; validation, S.M.-A. and R.G.; investigation, S.M.-A. and P.V.; visualization, S.M.-A.; writing—original draft preparation, all authors; writing—review and editing, all authors; supervision, P.G.-N.; project administration, P.G.-N. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by project PID2022-137334NB-I00 funded by MCIN/AEI/10.13039/501100011033 and by "ERDF/EU". This work has been partially funded by the Government of Aragón, through the research grant T32_23R Fluid Dynamics Technologies.

**Data Availability Statement:** The code and dataset for this work are fully available in Zenodo as a public repository: https://doi.org/10.5281/zenodo.13913312.

**Acknowledgments:** The authors would like to thank HYDRONIA LLC for the support provided for the development of this model using its software.
**Conflicts of Interest:** Authors S.M.-A. and P.G.-N., from University of Zaragoza, have a research agreement with the author R.G., from Hydronia HLL, for the development of new capabilities for the software commercialized by Hydronia HLL. The authors declare that the research was conducted in the absence of any commercial or financial interference that could be construed as a potential conflict of interest.

## References

1. International Tanker Owners Pollution Federation Statistics. 2024. Available online: https://www.itopf.org/knowledge-resources/data-statistics (accessed on 3 October 2024).

2. Chen, J.; Zhang, W.; Wan, Z.; Li, S.; Huang, T.; Fei, Y. Oil spills from global tankers: Status review and future governance. *J. Clean. Prod.* **2019**, *227*, 20–32. https://doi.org/10.1016/j.jclepro.2019.04.020. [CrossRef]

3. Spaulding, M.L. State of the art review and future directions in oil spill modeling. *Mar. Pollut. Bull.* **2017**, *115*, 7–19. https://doi.org/10.1016/j.marpolbul.2017.01.001. [CrossRef] [PubMed]

4. DHI. MIKE 21/3 Oil Spill: Oil Movement Modelling Software. 2024. Available online: https://www.dhigroup.com/technologies/mikepoweredbydhi/mike-21-3-oil-spill (accessed on 1 October 2024).

5. Shyue, S.w.; Sung, H.c.; Chiu, Y.f. Oil Spill Modeling Using 3D Cellular Automata For Coastal Waters. In Proceedings of the Seventeenth International Offshore and Polar Engineering Conference, Lisbon, Portugal, 1–6 July 2007.

6. Li, Y.; Lyu, X.; Ren, P. Oil Spill Timely Backtracking Oriented by Wind Field Correction With Self-Attentional Temporal Convolutional Networks. *IEEE J. Ocean. Eng.* **2024**, *49*, 114–132. [CrossRef]

7. Spaulding, M.; Kolluru, V.; Anderson, E.; Howlett, E. Application of three-dimensional oil spill model (WOSM/OILMAP) to Hindcast the Braer spill. *Spill Sci. & Technol. Bull.* **1994**, *1*, 23–35.

8. Zodiatis, G.; Lardner, R.; Alves, T.; Krestenitis, Y.; Perivoliotis, L.; Sofianos, S.; Spanoudaki, K. Oil spill forecasting (prediction). *J. Mar. Res.* **2017**, *75*, 923–953. [CrossRef]

9. Ivorra, B.; Gomez, S.; Ramos, A.M.; Glowinski, R. Nonlinear Advection–Diffusion–Reaction Phenomena Involved in the Evolution and Pumping of Oil in Open Sea: Modeling, Numerical Simulation and Validation Considering the Prestige and Oleg Naydenov Oil Spill Cases. *J. Sci. Comput.* **2017**, *70*, 1078–1104. [CrossRef]

10. Azevedo, A.; Oliveira, A.; Fortunato, A.; Bertin, X.; Bertin, A. Application of an Eulerian-Lagrangian oil spill modeling system to the Prestige accident: Trajectory analysis. *J. Coast. Res. J. Coast. Res. SI 2009*, *56*, 777–781.

11. Moghaddam, A.; Dabir, B. A 2-D hybrid particle tracking/Eulerian-lagrangian model for oil spill problems. *Indian J. Mar. Sci.* **2013**, *42*, 42–49.

12. Martínez-Aranda, S.; Ramos-Pérez, A.; García-Navarro, P. A 1D shallow-flow model for two-layer flows based on FORCE scheme with wet-dry treatment. *J. Hydroinform.* **2020**, *22*, 1015–1037. [CrossRef]

13. Murillo, J.; Martínez-Aranda, S.; Navas-Montilla, A.; García-Navarro, P. Adaptation of flux-based solvers to 2D two-layer shallow flows with variable density including numerical treatment of the loss of hyperbolicity and drying/wetting fronts. *J. Hydroinformatics* **2020**, *22*, 972–1014. [CrossRef]

14. Echeverribar, I.; Brufau, P.; García-Navarro, P. A fully Eulerian two-layer model for the simulation of oil spills spreading over coastal flows. *Phys. Fluids* **2023**, *35*, 116602. [CrossRef]

15. Martínez-Aranda, S.; Murillo, J.; Morales-Hernández, M.; García-Navarro, P. Novel discretization strategies for the 2D non-Newtonian resistance term in geophysical shallow flows. *Eng. Geol.* **2022**, *302*, 106625. https://doi.org/10.1016/j.enggeo.2022.106625. [CrossRef]

16. Ortega-Moya, J.; Martínez-Aranda, S.; Fernández-Pato, J.; García-Navarro, P. A vertically non-uniform temperature approach for the friction term computation in depth-averaged viscoplastic lava flows. *J. Comput. Phys.* **2024**, *519*, 113378. https://doi.org/10.1016/j.jcp.2024.113378. [CrossRef]

17. Echeverribar, I.; Martínez-Aranda, S.; Fernández-Pato, J.; García-Navarro, P. A GPU-based 2D viscous flow model with variable density and heat exchange. *Adv. Eng. Softw.* **2023**, *175*, 103340. https://doi.org/10.1016/j.advengsoft.2022.103340. [CrossRef]

18. Martínez-Aranda, S.; Murillo, J.; García-Navarro, P. A GPU-accelerated Efficient Simulation Tool (EST) for 2D variable-density mud/debris flows over non-uniform erodible beds. *Eng. Geol.* **2021**, *296*, 106462. [CrossRef]

19. Martínez-Aranda, S.; Fernández-Pato, J.; Echeverribar, I.; Navas-Montilla, A.; Morales-Hernández, M.; Brufau, P.; Murillo, J.; García-Navarro, P., Advances in Fluid Mechanics: Modelling and Simulations; Chapter Finite Volume models and Efficient Simulation Tools (EST) for Shallow Flows; Springer: Singapore, 2022; pp. 67–137.

20. AECOM Uganda Limited (AECOM). Tilenga Project Environmental and Social Impact Assessment (ESIA) Report. 2019. Available online: https://totalenergies.ug/tilenga-project-environmental-and-social-impact-assessment-esia-report (accessed on 1 September 2024).

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.