Modeling of the Fate and Behaviors of an Oil Spill in the Azemmour River Estuary in Morocco

Abstract

Oil spills are one of the most hazardous pollutants in marine environments with potentially devastating impacts on ecosystems, human health, and socio-economic sectors. Therefore, it is of the utmost importance to establish a prompt and efficient system for forecasting and monitoring such spills, in order to minimize their impacts. The present work focuses on the numerical simulation of the drift and spread of oil slicks in marine environments. The specific area of interest is the Azemmour estuary, located on Morocco’s Atlantic Coast. According to the environmental sensitivity index (ESI), given its geographical location at the intersection of the World’s Shipping Lines of oil transport, this area, as with many other sites in Morocco, has been classified as a high-risk area for oil spill accidents. By taking into account a range of factors, including the ocean currents, the weather conditions, and the oil properties, detailed numerical simulations were conducted, using the hydrodynamic TELEMAC-2D model, to predict the behavior and spread of an oil spill event in the aforementioned coastal region. The simulation results help to understand the spatial–temporal evolution of the spilled oil, the effect of wind on the spreading process, as well as the coastal areas that are most likely to be affected in the event of an oil spill accident. The simulations were performed with and without wind effects. The results showed that three days after the oil spill only 31% of the spilled oil remained on the sea surface. The wind was found to be the main factor responsible for oil drifting offshore. The results indicated that rapid action is needed to address the oil spill before it causes significant environmental damage and makes the oil cleanup process more challenging and expensive. The results of the present study are highly valuable for the management and prevention of environmental disasters in the Azemmour estuary area. The findings can be used to assess the efficacy of various response strategies, such as containment and cleanup measures, and to develop more effective emergency response plans.

1. Introduction

Estuaries are vital links between land and sea ecosystems, serving as important areas for leisure and economic activities. Estuaries are also ecosystems highly vulnerable to human-induced impacts and environmental changes, underscoring the need for adequate preservation [1]. The rise of maritime and river traffics, the high number of oil tankers, the construction of several submarine pipelines, as well as oil exploration in marine fields, have led to numerous oil spill accidents that pose a significant threat to estuaries and coasts [2]. Oil spills occur due to a variety of reasons, including oil exploration, ship collision, shipwrecking, oil unloading from tankers, ship tank cleaning, refinery activities, pipeline explosions, pipeline leaks, pipeline vandalism, and natural disasters [3,4]. Given its significant contributions to many government revenues, job creation, and energy supply, oil has been a crucial and highly sought-after resource in the world’s development [5]. According to a previous study by Chen et al. [6], the crude oil trade volume in 2017 reached 9947 billion ton-miles. Despite the relatively small proportion of oil that is lost in accidents, the total amount is very considerable. For example, tanker accidents alone caused a global loss of approximately 5.87 million tons of oil between 1970 and 2021 [7]. Oil spills are one of the most hazardous pollutants in marine environments.

The environmental impact of oil pollution has been extensively studied in numerous previous studies [8,9]. In their research on the Niger Delta community, Osuji and Onojake [8] discovered that hydrocarbons and heavy metals in crude oil have a significant negative impact on plants and animals by increasing the concentration of heavy metals in their cells. Ojimba et al. [9] investigated the socioeconomic repercussions of oil pollution on crop farms in Rivers State, Nigeria. They pointed out that crude oil pollution significantly reduced crop production and negatively affected the farmers’ income and their economic well-being. Overall, crude oil spills create a floating layer of plastic film in seawater, which hinders oxygen exchange with the air and prevents sunlight from entering the seawater, causing a significant decrease in plankton photosynthesis and growth. Since plankton are a vital component of the marine food chain, their decline has serious implications for the primary productivity of marine organisms [10]. Jabbar [11] examined how oil pollution affects the growth and diversity of aquatic plants in the Baher Al-Najaf Depression in Iraq. The author argued that oil pollution leads to the degradation of rare plant species and severely reduces plant distribution frequencies, affecting local food chain contamination. Allers et al. [12] found that coral colonies can be severely damaged or killed by sediments, including oil rig cuttings, despite their resilience to short-term sedimentation events. For example, the 2010 Gulf of Mexico oil spill (Macondo well) had a harmful impact on coral communities even in warmer deep waters. The oil spill also severely affected the abundance and composition of the bacterial communities in beach sands in the Gulf of Mexico. In addition, it affected the marsh vegetation of coastal saltmarshes in Barataria Bay, Louisiana, as well as the wetlands of the Mississippi River Delta, the farmed fish, and the shellfish operations in these wetlands. Moreover, it has had an impact on fish and shellfish farming in these wetlands [13].

When oil pollution occurs in icy areas, the layer of oil on the ice surface can absorb more solar energy, accelerating the ice-melting process [10]. However, it becomes difficult to locate oil spills in such regions because ice can hide them. Furthermore, the cold temperatures cause the oil to thicken and become more viscous (like tar), making it difficult to collect. According to Yang [14], cold environments are more vulnerable to oil spills because their ecosystems have adapted to harsh conditions, making them more susceptible to contamination. In cold environments, the most important factors in minimizing the environmental and economic impacts of oil spills are the ability and efficiency to contain and remove them, which vary from winter to summer. Various methods can be used to manage oil spills, such as physical/mechanical, chemical, biological, and sorbent treatments. Sorbent treatments are typically less expensive and commonly used to deal with oil spills. Recently, nanofiber sorbents have garnered attention from scientists as a potential material for absorbing oils [15,16].

Carpenter [13] asserted that international cooperation could help to reduce marine oil pollution. The authors also noticed that legal measures at regional, national, and international levels could increase the monitoring of discharges from ships. Moreover, the author indicated that surveillance activities using sophisticated available technologies, such as aerial and satellite surveillance, provide a good tool for environmental monitoring, enabling the identification of pollution sources and the prosecution of polluters. Through various surveillance and monitoring activities, measurement procedures, and technological development, oil pollution has considerably decreased in many parts of the world. However, underdeveloped countries and conflict zones may not have the same level of progress in reducing oil pollution.

Brakstad et al. [17] provided a concise summary of the behavior of crude oil spilled on the sea surface. As oil slicks drift, the volatile components of the crude oil evaporate into the air. Under the wave actions, part of the oil will be dispersed in the water column and water-in-oil emulsions will form on the sea surface. Over time, these emulsions will degrade and interact with suspended sediments, making them buoyant, particularly inshore or river areas with high-loading suspended sediment. This will be different with oil spills on the sea floor, such as the case of the Deepwater Horizon oil spill in the Gulf of Mexico at a depth of 1525 m, where the turbulence associated with the high-velocity discharge converted the oil, containing gas bubbles, into droplets with a wide range of sizes. The oil droplets will be carried upward through the water column, dispersing along a large sea surface. Understanding the behavior of oil spills is critical for modeling their impact.

Modeling is an essential tool for oil spill contingency planning, preparedness, and response. There are several approaches to modeling oil spills, often based on stochastic or deterministic methods. Over the past few decades, coastal and offshore oil spill research has expanded and evolved considerably. However, research on river oil spill has received less academic attention [18]. The stochastic approach involves a large number (tens to thousands) of individual deterministic simulations with varying input parameters, such as environmental conditions. The stochastic approach is typically used for risk assessment and preparedness planning. In contrast, the deterministic modeling approach simulates the release of oil from a single time point (a single set of wind and weather conditions). This method is commonly employed to predict the fate and behavior of oil using a single-operating model [19]. Niu et al. [19] examined the potential environmental effects of an oil spill in the Salish Sea located in western Canada. To better estimate the fate and behavior of the spilled oil in the target region, they employed a complete three-dimensional oil spill contingency and response (OSCAR) stochastic model. They also investigated the impact of chemical dispersant use. The results of this study have the potential to enhance the accuracy and efficiency of oil spill response planning, including decision support on chemical dispersant application, as well as aid in forecasting environmental impacts.

Current entrainment and diffusion of a tracer, considering source and sink terms, such as oil spills, can be also simulated using the TELEMAC model. This model can predict the quality of seawater and the dispersion of a plume pollutant in 1D, 2D, and 3D, taking into account the weather conditions and local currents in a particular area. Several studies were conducted using the TELEMAC model to simulate various hydrodynamic ocean conditions with and without oil spills [20,21,22,23,24,25]. Lopes et al. [20] employed the TELEMAC-3D hydrodynamic model coupled with the Easy Coupling Oil System (ECOS) model to simulate oil spill behaviors under different extreme environmental conditions in the Patos Lagoon, located in the southern region of Brazil. The TElEMAC-3D model proved capable of replicating the current patterns in the estuary of Patos Lagoon, leading to the development of effective oil spill response plans under flood and ebb conditions in the study area. To investigate the impact of sea level rise (SLR) on oil spill dispersion, Lavine et al. [21] conducted a numerical study using the hydrodynamic TELEMAC-2D model to predict oil spill spreading in the Pulai River estuary and southwest Johor Strait before and after the SLR phenomenon. The simulated results demonstrated good agreement with the measured values obtained. Eke et al. [22] argued that oil spills in estuaries have received less attention and are not as fully understood as their oceanic counterparts. In an attempt to address this gap, the authors provided a detailed analysis of estuarine oil spill transport for the Humber Estuary in the UK, using a coupled hydrodynamic TELEMAC-3D model and oil spill model. The research study produced important and novel findings on how the interaction of river discharge and tidal range influences oil slick dynamics, which could aid responders in assessing likely oil trajectories and deploying response tools in a cost-effective manner. Monteiro et al. [23] noted that the variability of the oceanographic and meteorological conditions are sources of uncertainties in the susceptibility prediction of oil contamination. The authors suggested that the reliability of susceptibility assessments could be improved by incorporating a greater number of oil spill scenarios into the analysis.

Zhang [26] reported that operational oil spill monitoring is currently involved using a combination of optical and radar satellite imageries. The author noted that the Global Positioning System (GPS), Global Navigation Satellite System (GNSS), Galileo system, and BeiDou system provide good signal sources to monitor the ocean. To detect oil-slicked ocean surfaces, Zhang [26] used coastal simulation based on a real oil spill accident that caused marine oil pollution, using BeiDou Inclined Geosynchronous Orbit (IGSO) and Medium Earth Orbit (MEO) MEO satellites. The simulation involved calculating a delay-Doppler map (DDM) of the oil-slicked surface using a mean-square slope (MSS) model for oil-slicked/clean surfaces and a Zavorotny–Voronovich (Z-V) scattering model. The results showed that the oil-slicked area can be clearly distinguished by using DDM map technology from BeiDou MEO and IGSO satellites.

Kvočka et al. [18], reviewing the state of the art of river oil spill modeling and summarizing developments in the field from 1994 to 2021, revealed that many gaps in knowledge still exist. The authors suggested that there is a need for: (i) experimental studies to calibrate and validate numerical models and better understand the main physicochemical processes, (ii) the governing processes, i.e., hydrodynamics, advection–dispersion, weathering processes, should be interlinked, (iii) adaptation and validation of coastal and offshore oil spill models for fluvial environments, and (iv) development of remote sensing systems and techniques for detecting oil spills in rivers. The authors highlighted the difference in the modeling approach for river and offshore oil spills due to many factors, i.e., water density, flow hydrodynamic structures (rapid change of flow regimes in rivers), velocity distributions (dominance of downstream current in rivers), presence of vegetation and more interaction with spilled oil in river environments, more frequent interaction of oil with sediments in river systems, and turbulent mixing overall flow depth of a river. During the discussion about simulating oil spills in river environments, different numerical models with one, two, and three dimensions were reviewed. Examples of these models include TELEMAC-2D (a 2D depth-averaged model), CE-QUAL-W2 (a 2D water quality and hydrodynamic model used in stratified surface water systems in reservoir, lake, and estuarine systems), GeoSpatial Stream Flow Model (GeoSFM), Subsurface Oil Simulator (SOSim) model, HEC-RAS (a 1D hydraulic model), ADH (a 2D ADaptive Hydraulics model), EFDC (a 3D Eulerian hydrodynamic model), and GLLVHT (a 3D Generalized, Longitudinal–Lateral–Vertical Hydrodynamic, and Transport model). Most of these models have the capability to be coupled with oil/chemical spill models and can be dynamically linked in real-time with streamflow and meteorological data.

This study aims to understand the potential environmental damage and identify the marine and coastal areas most likely to be affected in the event of an oil spill accident in the Azemmour estuary. To achieve this, the TELEMAC model was used to simulate the dispersion and diffusion of an oil slick, and the outcomes were utilized to design more accurate and efficient plans for responding to oil spills. These plans will include strategies for containing and cleaning up the oil spill, as well as identifying the areas most at risk of environmental damage. The use of simulated outcomes can help mitigate the environmental impact of oil spills and reduce harm to marine and coastal ecosystems. The ultimate goal of this study is to safeguard the ecosystem of the Azemmour estuary and minimize the consequences of any potential oil spill accidents.

2. Study Area

This study focuses on the primary coastal area situated at the estuary of the Oum Er-Rbia River, as shown in Figure 1. This area is located on Morocco’s Atlantic Coast in the Azemmour region, around 80 km southwest of Casablanca. The coastal area is characterized by the presence of sand dunes about 1000 m wide, mostly formed by north–northeast winds that are typical of the Moroccan coastal zone. Both El Haouzia and Lalla Aicha El Bahria beaches, located near the river mouth, have pebble deposits that are a mixture of sedimentary and eruptive origin [27]. The river mouth in the estuary has a variable width of approximately 150 m.

The study area is characterized by steep shores and tree-covered cliffs. The river section that connects the ocean extends straight for about 3 km, and its banks are flat, bare, or covered with shrubs, creating an environment conducive to the development of tidal bores upstream of the river mouth. These tidal bores can extend over a significant distance, and they have a considerable impact on the estuary’s ecosystem.

The hydrodynamic structures in the estuarine area are primarily influenced by the river flow and tidal currents. The construction of several dams along the Oum Er-Rbia River has led to a significant reduction in freshwater inflows. Usually, the highest monthly mean flow of the Oum Er-Rbia River occurs in March, while the lowest one occurs in August with a flow rate of almost 40 m³/s. Along the Atlantic Moroccan coast, the tide is semi-diurnal, and at the estuary of the Oum Er-Rbia River, its amplitude alternates between approximately 1.5 m and 4 m, as reported by the Casablanca Port maritime services [27].

3. Model Theoretical Formulation

The hydrodynamic model employed in this study is the TELEMAC-2D model. This model solves the depth-averaged free-surface flow equations of Saint-Venant based on the finite-element or finite-volume methods and a computation mesh of triangular elements. It is able to simulate both transient and permanent conditions [28]. The model accounts for various phenomena, such as the propagation of long waves, bed friction, Coriolis effect, meteorological parameters, flow properties and turbulence, water density variation, entrainment, and diffusion of a tracer by currents, drag forces generated by vertical structures, and many other parameters. The TELEMAC-2D model simultaneously solves the following four hydrodynamic equations:

Continuity equation

$$\frac{\partial d}{\partial t}+V\nabla \left(d\right)+ddiv\left(\overrightarrow{V}\right)=S_d$$

(1)

Momentum equation along the x-direction

$$\frac{\partial u}{\partial t}+V\nabla \left(u\right)=-g\frac{\partial \eta }{\partial x}+S_X+\frac{1}{d}div\left(d{\mathsf{\nu }}_t\nabla \left(u\right)\right)$$

(2)

Momentum equation along the y-direction

$$\frac{\partial v}{\partial t}+V\nabla \left(v\right)=-g\frac{\partial \eta }{\partial y}+S_y+\frac{1}{d}div\left({d\mathsf{\nu }}_t\nabla \left(u\right)\right)$$

(3)

Tracer conservation equation

$$\frac{\partial T_p}{\partial t}+V\nabla \left(T_p\right)=S_T+\frac{1}{d}div({d\mathsf{\nu }}_T\nabla \left(T_p\right))$$

(4)

where d is the water depth, t is the time, V = (u, v) is the depth-averaged flow velocity vector, u and v are the depth-averaged velocity components in the x- and y-directions, (x, y) are horizontal space coordinates, ∇ is the gradient operator, div is the divergence operator, S_d is a source or sink of fluid, g is the gravitational acceleration, η is the free surface elevation, (S_x, S_y) are the source or sink terms in the dynamic equations representing the wind, S_T is the source or sink of tracer, (ν_t, ν_T) are momentum (eddy viscosity) and tracer diffusion coefficients, and T_p is the passive (non-buoyant) tracer. Herein the unknowns are d, u, v, and T_p.

In this study, the computational domain size is not large enough, allowing for the generation of a high-resolution mesh, and therefore, a k-Epsilon model was used to describe the turbulent flow motions in the horizontal (x, y)-plane. The depth-averaged turbulent kinetic energy (k) and its relative dissipation rate (ε) are expressed by the following equations:

$$\frac{\partial k}{\partial t}+V\nabla \left(k\right)=\frac{1}{d}\nabla \left(d\frac{{\mathsf{\nu }}_t}{{\sigma }_k}\nabla \left(k\right)\right)+P-\epsilon +P_{kv}$$

(5)

$$\frac{\partial \epsilon }{\partial t}+V\nabla \left(k\epsilon \right)=\frac{1}{d}\nabla \left(d\frac{{\mathsf{\nu }}_t}{{\mathsf{\sigma }}_{\epsilon }}\nabla \left(\epsilon \right)\right)+\frac{\epsilon }{k}\left(C_{1\epsilon }P+C_{2\epsilon }\epsilon \right)+P_{\epsilon v}$$

(6)

with

$$P={\mathsf{\nu }}_t\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\frac{\partial u}{\partial y}$$

(7)

$$P_{kv}=C_k\frac{u_{*}^3}{d}\hspace{1em}\mathrm{a}\mathrm{n}\mathrm{d}\hspace{1em}{P}_{\epsilon v}=C_{\epsilon }\frac{u_{*}^4}{d^2}$$

(8)

where C_k = C_f^(−1/2) and C_ε = 3.6(C_2εC_µ^(1/2))/C_f^(−1/2), C_f is the bed friction coefficient, u_∗ is the shear velocity, ν_t = C_µ(k²/ε), C_µ = 0.09, C_1ε = 1.44, C_2ε = 1.92, σ_k = 1.0, σ_ε = 1.3. The use of the k-Epsilon model often requires a finer mesh than the constant viscosity model and, in this way, increases computation time. The model equations are solved by a fractional step method, with the convection of turbulent variables being processed at the same time as the hydrodynamic variables, and the other terms relating to the diffusion and production/dissipation of turbulent values being processed in a single step.

The hydrodynamic TELEMAC-2D model will then be coupled with a transport system model to simulate the oil spill behaviors under various environmental conditions. The transport model simulates different processes, including the oil slick advection and diffusion as well as the effect of weather conditions such as the oil spreading, evaporation, dissolution, and volatilization. Goeury et al. [29] proposed detailed theoretical approaches for these processes. Here, for simplicity, we briefly mention some of the oil spill mechanisms.

The drifting of an oil slick on a free water surface is primarily induced by both the current and wind velocities. Thus, the drift velocity of the oil slick can be expressed as:

$$V_o=V_c+\mathsf{\beta }{V}_w$$

(9)

where V_o is the oil-slick velocity vector, V_c is the current velocity vector at the free surface, V_w is the wind velocity vector, and β is a wind drift factor. The current velocity can be estimated using the classical log-law assumption based on the depth-averaged velocity. By applying the second law of dynamics to the advected oil particle in steady condition, and after some assumptions and simplifications, considering the relative particle velocity with respect to the current and wind velocities, a value of almost 3.6% was defined for β.

The oil particle trajectory is considerably affected by flow turbulence. To account for this, a stochastic approach was adopted based on a bi-dimensional advection–diffusion equation transformed into a Lagrangian representation, Equation (4), interpreting the oil concentration as a probability of finding a particle at a given (x, y)-location and time t. Goeury et al. [29] proposed a stochastic solution to the abovementioned equation as follows:

$$\mathrm{X}\left(t+\delta t\right)=X\left(t\right)+\left[V_o-\frac{1}{d}\nabla \left(\frac{d{\nu }_t}{{\sigma }_c}\right)\right]\delta t+\sqrt{\frac{2{\nu }_t}{{\sigma }_c}\delta t}\mathsf{\xi }\left(t\right)$$

(10)

where X(x, y) is the particle location, σ_c is the neutral turbulent Schmidt number and ξ(t) is a vector with independent, standardized random components.

Typically, an oil spill on a water surface spreads quickly over a large area, where the spreading process is often controlled by the gravitational, inertial, surface tension and viscous forces. As most of the surface of an oil spill (>90%) is controlled by gravity, the spreading process can be simply predicted as follows [29]:

$$S={\left(\frac{27\pi }{2}\frac{W^{3}g∆}{{\nu }_o}t\right)}^{\frac{1}{4}}$$

(11)

where S is the slick surface, W is the volume of spilled oil, Δ = (ρ_w − ρ_o)/ρ_w, ρ_w is the water density, ρ_o is the oil density, and ν_o is the oil kinematic viscosity.

An oil spill in the ocean loses a significant amount of its initial volume over time by evaporation, which represents the most crucial mass transfer process in such an oil spill. The change in mass of an oil component i is modelled [29] as:

$$\frac{d{m}_i}{dt}=-K_{ev}{A}_i\frac{m_{i}exp\left[\frac{∆H_i}{R{T}_{Bi}}\left(1-\frac{T_{Bi}}{T}\right)\right]}{\sum _j\frac{m_j}{M_j}RT}$$

(12)

where m_i is the mass of component i, K_ev is the evaporation mass transfer coefficient, A_i is the area of component i, ΔH_i is the molar enthalpy of component i, R is the universal gas constant, T_Bi is the boiling point of component i, T is the ambient temperature, m_j is the mass of an element j, and M_j its molar mass.

Other algorithms to estimate the dissolution and volatilization processes were also proposed. The characteristic coefficients of the different processes involved in the transport model were determined through laboratory experiments. Verification of the numerical model was investigated using benchmark test cases and well-recorded real oil spill cases.

4. Model Boundaries and Forcing Data

4.1. Bathymetry

The bathymetry of the study area was generated using the SHOM charts, SHOM6120 (source: https://www.nauticalchartsonline.com/charts/SHOM/Morocco, accessed on 3 May 2023), Lat. 33°23′, with a scale of 1:50,027. Figure 2 displays an enlarged section of the Bathymetry of the coarse grid domain at the Azemmour estuary.

Figure 2 shows that the coastal zone of the study area has relatively regular isobaths that run parallel to the shoreline. From the shoreline to the −10 m bathymetric line, a smooth bottom slope of 0.5% extends for about 2.2 km. Between the bathymetric lines of −10 m and −20 m, the slope of the bottom softens to an average value of 0.1%. The bathymetric line at −20 m is nearly 9 km from the coastline.

An accurate representation of boundary geometry is a crucial component for successful numerical modeling. The bathymetry was extracted in (x, y, z)-points and a computational mesh was generated for the spatial domain. The Blue Kenue software was used to generate unstructured and regular triangular meshes.

The resulting computational domain consists of an unstructured triangular mesh with an edge length of 50 m, covering an x- and y-length of approximately 8 km and 7 km, respectively (see Figure 3). Previous research [30,31] suggested that a high-resolution mesh, ranging from 20 to 100 m edge length, captures the detailed bathymetric feature, which is sufficient to model the flow motion and wave propagation accurately. Additionally, a high-resolution mesh is not expected to significantly affect the numerical results, as reported by Eke et al. [22].

4.2. Meteorological Parameters

Accurate weather data is critical to mapping oil spill vulnerability and conducting accurate numerical simulations [3]. Environmental factors such as winds, waves, and currents play a significant role in determining the movement of an oil spill. In this study, the relevant meteorological parameters were analyzed. The meteorological data were collected from the weather station closest to the Azemmour estuary, located at El Jadida, 15 km southwest of the study site. Figure 4 illustrates the annual average wind direction rose. Figure 4 indicates that the most frequent wind direction in the area is from the north–northeastern sector, with an annual average velocity of less than 5 m/s. Due to limited spatial and temporal wind data in the study area, numerical simulations were conducted based on an average wind velocity of almost 4 m/s over the computational domain.

In the study area, the hydrodynamic model accurately reproduced the depth-averaged currents by incorporating precise bathymetry and average tidal height data at the offshore boundary. The tide propagates from the southwest and has a semi-diurnal period of 12 h 25 min and an average amplitude of 2.8 m, as shown in

Table 1. Tidal levels in the Azemmour estuary.

Tide	Low Tide	High Tide	Tidal Range
Exceptional high water	0.20	3.90	3.70
Medium–high water	0.70	3.50	2.80
Medium still waters	1.50	2.80	1.30

.

5. Results

5.1. Hydrodynamic Model Results

In order to conduct oil spill modeling, it is necessary to have an accurate understanding of the hydrodynamic conditions of the sea, which can be obtained from surveys, measurements, or simulations based on characteristic parameters derived from historical data of the target area. However, conducting surveys and measurements over long periods of time to describe sea conditions can be expensive and not always feasible. As an alternative, numerical modeling can be used to describe the sea’s hydrodynamic conditions, including surface elevation and sea current components, over large scales and long periods. This approach can provide a cost-effective and practical solution to obtaining the necessary data for oil spill modeling [32].

In the present study, the hydrodynamic model was developed, taking into account various driving forces, for annual average weather conditions in the Azemmour estuary. The model considers a typical annual mean flow discharge of 117 m³/s at the mouth of the Oum Er-Rbia River. Since the computational domain size is not large enough, the simulations were performed under barotropic conditions, with constant water temperature and salinity. At the offshore open boundaries, a tidal height of almost 2.8 m (see Table 1) was imposed, and a solid condition is used for the shorelines. Due to the unavailability of spatial wind data in the study area, a wind velocity of constant direction and magnitude of almost 4 m/s was considered. The numerical simulation was performed over a period of three days, which was considered sufficient to provide useful insights into the behavior of the oil spill in such a computational domain, as shown in Figure 3.

Figure 5 shows the mean water depth, d, in the computational domain of the study area, averaged over the simulation period. The water depth gradually increases from the mouth of the Oum Er-Rbia River going offshore in relatively regular iso-depth bands parallel to the shoreline, which corresponds with the bathymetric distribution. This regular distribution of water depth could also allow for regular sea current distribution.

Figure 6 illustrates the distribution of the sea current speed U in the x-direction across the computational domain. The plot shows that U remains nearly constant over significant portions of the computational domain, with values mostly ranging between 0.6 and 0.8 cm/s, reaching a maximum of almost 1 cm/s near the coastline. Figure 6 clearly shows a smooth and gradual variation of U in the x-direction, reflecting the effect of the gradual flow depth, as shown in Figure 5. In Figure 6, positive values of U are observed across most of the computational domain, while a local vortical structure development is shown in the upper right corner of the plot, where U takes on negative values (indicated by the red color).

Figure 7 shows the distribution of the computed current speed V in the y-direction. Similar to the U-component, V also exhibits positive values in almost all computational areas. Notably, the current magnitude in the y-direction is almost seven times greater than that in the x-direction, with speeds ranging from 4.5 to 5.5 cm/s. This indicates that the sea current in the study area is primarily directed toward the north-northeast, as illustrated by the velocity vector map in Figure 8. Additionally, in the upper right corner of Figure 7, V also shows small or negative values, which confirms the development of a local vortical structure in this location.

The predicted U and V velocities of the sea current, obtained from the hydrodynamic model, will be used later as input to predict the spreading of the spilled oil (see Section 3).

5.2. Oil Spill Results

To model an oil spill in the marine environment, it is necessary to calibrate the dissolution and volatilization mass transfer coefficients, which can be achieved through artificial tests that take into account the properties of the hydrocarbon and the mean current velocities. The resulting hydrocarbon kinetics were then used to calibrate the numerical transport model, as described in [29].

It is important to note that the outcome of an oil spill can be influenced by various factors, such as the type and quantity of oil spilled, the distance from the release location, and weather conditions. In this study, the mass of the oil spill considered at the beginning of the simulation was approximately 1700 kg. The location of the oil release was selected near the Oum Er-Rbia River’s mouth, at coordinate (x = 559,000, y = 3,688,000), indicated by the star symbol in the figures below. This location was chosen on the shipping lane to the port of Al Jadida, subject to a high risk of oil accidents, and is located in a site of high ecological, environmental, and economical values. The simulated oil type in this study is moderate crude oil.

The numerical simulation revealed that after three days, only around 31% of the initial oil spill volume remained on the ocean surface, far from the shoreline. The amount of oil stranded on the shoreline, adhering to substrate surfaces, was about 57%, while roughly 12% of the oil mass had evaporated. The amount of dissolved oil mass in the water column was relatively insignificant, accounting for only about 0.1% of the spilled oil mass. These results are attributed to the relatively calm weathering conditions in the Azemmour estuary.

Oil spilled on the water surface can spread in all directions, but the sea current, flow turbulence, and winds can cause the oil slick to advect in a particular direction. Figure 9 shows that, after three days, most of the remaining oil on the surface has moved further northwest, while a portion is moving southwestwards. Although the simulated sea currents are north northeastwards, as shown in Figure 8, the wind forces acting on the water surface strongly influence the oil slick’s movement. It is worth mentioning that the simulation used a northeast wind of greater magnitude than the current speeds. The oil trajectory, illustrated in Figure 9, shows a northeast direction, reflecting the effect of the sea current on the oil dispersion. These results indicated that the spilled oil can quickly spread over a large water surface and travel a significant distance offshore. Prompt action is necessary to address the situation before it causes significant environmental damage. As time passes, the cleanup process becomes more challenging and expensive. Given that around 57% of the oil can remain stranded, an immediate response strategy is crucial for cleaning up spilled oil. Oil reaching the shoreline can cause substantial damage to the local fauna and flora, coastal communities, and tourist infrastructure.

Figure 10 shows the results of the numerical simulation of the spilled oil spread, without considering the wind effect. The simulation shows that, after three days, the oil has spread over a vast area, from the shoreline to the offshore side. The ocean current, with its main northeasterly direction, is the only factor influencing the simulated oil movement. The results also show that most of the remaining oil on the water surface is moving towards the shore areas. This could result in a large mass of oil reaching the shorelines and flowing upstream along the mouth of the Oum Er-Rbia River, causing considerable damage to the estuary ecosystem. These findings emphasize the critical need for immediate action to prevent the further spread of oil, saving lives and property.

6. Conclusions

Oil spills can cause extensive damage to marine ecosystems and biodiversity. The spread and dispersion of oil spills are mainly influenced by seawater currents and winds. Understanding these processes is crucial for effective oil spill cleanup and response planning, as well as predicting environmental impacts. Fortunately, recent advances in software, such as TELEMAC, have significantly improved the ability to model and predict oil spill motion. This powerful modeling tool allows engineers and scientists to simulate the fate and behavior of an oil spill under various weather and oceanic conditions.

The purpose of this study was the numerical simulations of an oil spill scenario in the Azemmour estuary in Morocco, using the TELEMAC-2D model. The simulations took into account the impact of various environmental conditions on the fate and behavior of the spilled oil.

The study was conducted in two phases. In the first phase, a two-dimensional hydrodynamic modeling was conducted to understand the circulation of marine currents and the flow turbulence levels in the study area. The results of this modeling were subsequently used in an oil spill model to predict the behavior of the spilled oil.

The results of the numerical simulation revealed that, after three days of the oil spill, only 31% of the spilled oil remained on the sea surface. The wind was found to be the main factor responsible for drifting the oil offshore, while the ocean current’s direction influenced its trajectory. The rapid spread of spilled oil highlights the need for immediate action to prevent severe environmental and economic damage that oil spills can cause.

In the absence of wind effect, the simulations revealed that, under the current driving forces, a large amount of spilled oil could reach the shorelines and move upstream along the Oum Er-Rbia River’s mouth, leading to severe ecological damage and potential economic and health consequences. Therefore, immediate action to prevent the oil spill is crucial.

Numerical simulations, such as those performed in this study using TELEMAC, are essential for predicting sea current circulations and their effects on spilled oil dispersion. The findings of the present study can be used to evaluate the effectiveness of different response strategies, develop emergency response plans and guide future research, such as modeling more specific weather conditions and installing coastal monitoring stations for more accurate measurements.

For the Azemmour estuary, it is of paramount importance to continue research and simulations on oil spills to better understand the effects of different environmental conditions and to improve response strategies. Further simulation scenarios should include more specific weather conditions such as seasonal variation, different oil properties, different locations of oil spilling, and variable wind distribution. In addition, installing coastal monitoring stations and carrying out field current measurements are highly recommended to provide more accurate and real-time data on ocean currents, improving the accuracy of modeling and prediction of oil spill behavior. A more detailed analysis of the dynamics of oil spills and contaminant migration throughout the study area is also highly suggested. Furthermore, a 3D modeling approach is recommended to obtain a more complete description of the hydrodynamic process, which can aid in the development of more effective emergency response plans and mitigation efforts.