Abstract

Atmospheric dispersion modelling in West Africa is crucial to understand, predict, and mitigate the impacts of air pollution on public health, ecosystems, and regional climate in a context of rapid urbanization and industrialization. This study develops a three-dimensional analytical model for carbon monoxide (CO) dispersion, governed by the turbulent scalar flux transport equation within a second-order Eulerian framework. A theoretical analysis of the model establishes the existence and uniqueness of the solution, and derives the associated optimality system using control theory and monotone operator techniques. From this formulation, an analytical solution is obtained that explicitly incorporates diffusion coefficients, turbulence intensities dependent on atmospheric stability (cloud cover and radiation), and three-dimensional wind components $\left(u,v,w\right)$ derived from ERA5 reanalysis. The analytical model is then calibrated through four-dimensional variational data assimilation (4D-Var) using satellite-based CO observations from the Copernicus EAC4 reanalysis, to produce fields consistent with observations. An application over the Gulf of Guinea during the dry season demonstrates that the coupled Model-4D-Var system generates CO dispersion maps with finer spatial resolution (0.25˚ × 0.25˚) than the input satellite products (0.75˚ × 0.75˚), while maintaining good agreement with observations, as indicated by low RMSE, small negative fractional bias, and FAC2 values between 0.5 and 2 at most pressure levels. These high-resolution maps enable the quantification, identification, and tracking of CO plumes in space and time, including persistent plumes over southern and northern Nigeria, Ghana, and the Central African Republic, linked to gas flaring, biomass burning, and urban/industrial emissions. The results highlight the potential of combining an analytical dispersion model with 4D-Var to improve CO monitoring in data-sparse regions, such as the Gulf of Guinea.

Keywords

Carbon Monoxide, 4D-Var Data Assimilation, ERA5, EAC4, Gulf of Guinea

Share and Cite:

1. Introduction

Air pollution is a leading cause of morbidity and premature mortality worldwide and represents a major environmental and public health challenge at the global scale [1]. Numerous anthropogenic and natural sources contribute to atmospheric pollution, including fossil fuel combustion, biomass burning, gas flaring, forest fires, road traffic, residential heating, and industrial activities [1]-[4]. Emissions from biomass burning, in particular, contain large amounts of gaseous compounds and particles that significantly affect local, regional and even global air quality [5]. Among the emitted species, carbon monoxide (CO) plays a key role in controlling the oxidative capacity of the troposphere through its reaction with the hydroxyl radical (OH), and thereby influences ozone formation and loss depending on nitrogen oxide (NO_x) levels [6]. Typical CO lifetimes of several weeks allow long-range transport and regional accumulation, especially in the tropics [7].

In West Africa, and more specifically over the Gulf of Guinea, CO originates from both natural and anthropogenic sources, including agricultural fires, savanna and forest fires, domestic combustion, and gas flaring associated with oil and gas production [4] [8]-[10]. Nigeria is among the most affected countries due to intensive flaring and rapid urbanization, with cities such as Onitsha and Kaduna reported among the most polluted worldwide in terms of particulate matter [1] [11] [12]. Other coastal and inland countries (Ghana, Côte d’Ivoire, Burkina Faso, Niger, Togo, Benin) also contribute through industrial emissions, road transport, domestic fires and biomass burning [10] [13]. In addition, long-range transport of dust from the Sahara-Sahel region and biomass-burning plumes from Central Africa further enhances pollution levels over the Gulf of Guinea, particularly during the dry season (December-February) [8] [13] [14].

Despite this strong pollution burden, air quality modeling over West Africa remains among the least developed in the world, owing to limited observational networks, scarce emission inventories, and reduced computational and financial resources [12] [13]. Atmospheric dispersion modeling is inherently complex, as it requires accounting for a wide range of physical and chemical processes, including advection, turbulent diffusion, boundary-layer dynamics, and, where relevant, chemical transformation [15]. While several observation networks and campaigns (e.g. AERONET, AMMA, DACCIWA) have improved our understanding of aerosols and trace gases in the region [8] [13], few studies have addressed the explicit modeling of CO dispersion at regional scale over the Gulf of Guinea [4] [16]. Understanding the spatio-temporal dynamics of CO plumes in this region, therefore remains a challenging and partly open problem.

In this study, we focus on quantifying and characterizing CO dispersion over the Gulf of Guinea using a physically-based yet computationally efficient modeling framework. We develop a three-dimensional analytical dispersion model derived from the complete turbulent scalar flux transport equation (Differential Flux Model, DFM) [17], which accounts for horizontal and vertical diffusion under stability-dependent turbulence parameterizations. The model is driven by ERA5 reanalysis meteorological fields [18] and is coupled with a four-dimensional variational (4D-Var) data assimilation scheme used to adjust the initial CO concentration based on satellite-derived CO data from the ECMWF Atmospheric Composition Reanalysis 4 (EAC4) [19].

The main contributions of this work are threefold. First, we derive and analyze a 3D analytical solution for CO concentration, including a theoretical study of existence, uniqueness and optimality using monotone operator theory and optimal control [20]-[22]. Second, we implement a 4D-Var assimilation framework that constrains the analytical model with ERA5 meteorology and EAC4 CO, and we evaluate its performance using standard statistical metrics (RMSE, fractional bias, FAC2) [23] [24]. Third, we apply the coupled system over the Gulf of Guinea to identify and track major CO plumes associated with gas flaring, urban emissions, and biomass burning, and we discuss the role of regional wind regimes and atmospheric stability in shaping their dispersion [2] [4] [8].

The paper is organized as follows. Section 2 describes the study area, data and methods, including the analytical formulation and the 4D-Var setup. Section 3 presents the numerical experiments, validation results, and analysis of CO plumes. Section 4 summarizes the main findings and outlines perspectives for future improvements and applications.

2. Study Area, Data and Methods

2.1. Study Area

The study focuses on coastal areas, specifically in the Gulf of Guinea (Figure 1), which extends from 4˚S to 18˚N and from −10˚W to 20˚E. It represents 50% of the continent’s pollution and 10% globally [11]. The study area is characterized by two seasons: the rainy season and the dry season. The long dry season extends from December to February; the short dry season, from July to September. As for the rainy season, the long rainy season occurs from May to June and begins in March-April; the short rainy season, from October to November [25]. The West African region is characterized by four types of wind regimes: the monsoon and Harmattan winds, the nocturnal low-level jet (NLLJ) and the land-sea breeze (LSB) system [25]. These different wind regimes play a crucial role in the redistribution of pollutants emitted by coastal area cities [25].

Figure 1. The image shown here delimits the study area by highlighting coastal countries with a black line.

However, we denote $\mathcal{Q}$ the study area with a spatial resolution of 0.25˚ × 0.25˚ (approximately 28 km × 28 km) and a temporal resolution of three hours. It is distributed according to latitude, longitude, and pressure level (altitude). Each point on the grid, according to the pressure level, has an area of 28 km × 28 km.

2.2. Data

In this study, we used meteorological data from the ERA5 reanalysis developed by the European Centre for Medium-Range Weather Forecasts (ECMWF) [18]. This data is accessible via the following link: https://cds.climate.copernicus.eu/. It is characterized by a spatial resolution of 0.25˚ × 0.25˚ and three-hourly temporal resolution. These meteorological data are presented in Table 1. TOA incident radiation, fraction of cloud cover and friction velocity are defined without pressure levels.

Table 1. Meteorological variables (ERA5).

As for the carbon monoxide (CO) data, they come from the EAC4 (ECMWF Atmospheric Composition Reanalysis 4; https://atmosphere.copernicus.eu) [19]. They are extracted at a resolution of 0.75˚ × 0.75˚, at a time resolution of 3 hours and over several atmospheric pressure levels.

2.3. Methods

The modeling of pollutant dispersion in our area is carried out through the combined application of two methods. The first methodological approach consists of developing a mathematical model that represents the dynamics of pollutant dispersion. The second step involves calibrating the analytical model, which is derived from the mathematical model.

2.3.1. Mathematical Formulation

In the context of this study, we consider $\mathbb{Q}$ as the theoretical domain of study of the model such that:

$\{\begin{array}{l}\Omega =\left\{ℝ,ℝ,{ℝ}^{+}\right\},\\ \left(x,y,z,t\right)\in \mathbb{Q}=\Omega \times \left[0,T\right[,\\ \left(x,y,z,0\right)\in {\mathbb{Q}}_0=\Omega \times \left\{0\right\},\\ \Lambda =\partial \Omega \times \left[0,T\right[.\end{array}$ (1)

The dynamics of pollutant dispersion are governed by the equation:

$\begin{array}{l}\frac{\partial }{\partial t}\left(K_x\frac{\partial c}{\partial x}\right)+\frac{\partial }{\partial t}\left(K_y\frac{\partial c}{\partial y}\right)+\frac{\partial }{\partial t}\left(K_z\frac{\partial c}{\partial z}\right)-V_x\frac{\partial }{\partial x}\left(K_x\frac{\partial c}{\partial x}\right)\\ -V_y\frac{\partial }{\partial y}\left(K_y\frac{\partial c}{\partial y}\right)-V_z\frac{\partial }{\partial z}\left(K_z\frac{\partial c}{\partial z}\right)\\ =-\left({\sigma }_u^2\frac{\partial c}{\partial x}\right)-\left({\sigma }_v^2\frac{\partial c}{\partial y}\right)-\left({\sigma }_w^2\frac{\partial c}{\partial z}\right)+\frac{\partial }{\partial x}\left[{\sigma }_u{\wedge }_x\frac{\partial }{\partial x}\left(K_x\frac{\partial c}{\partial x}\right)\right]\\ +\frac{\partial }{\partial y}\left[{\sigma }_v{\wedge }_y\frac{\partial }{\partial y}\left(K_y\frac{\partial c}{\partial y}\right)\right]+\frac{\partial }{\partial z}\left[{\sigma }_w{\wedge }_z\frac{\partial }{\partial z}\left(K_z\frac{\partial c}{\partial z}\right)\right],\end{array}$ (2)

where $c$ is the pollutant concentration; $V_x$ , $V_y$ and $V_z$ are the components of the wind speed vector $u$ ; ${\sigma }_u$ , ${\sigma }_v$ and ${\sigma }_w$ are the standard deviations (or the variances of the along-wind and horizontal and vertical accross-wind speed respectively) of turbulent velocity fluctuation in the $x$ , $y$ and $z$ directions, respectively. $K_x$ , $K_y$ and $K_z$ are a diffusivity coefficient (in units of m²/sec); ${\wedge }_x$ , ${\wedge }_y$ and ${\wedge }_z$ corresponds to the diffusion scale which is a macroscopic scale [26]. The main Equation (2) is derived from the transport equation for turbulent scalar fluxes [17]. This equation includes first- and second-order tensors. Closing the first-order tensors yields the diffusion coefficients [27]. Similarly, closing the second-order tensors yields the velocity standard deviations and diffusion scales [27]. We assume the following assumptions:

1) The diffusion coefficients, velocity coefficients, and diffusion scales are constant;

2) At $t_0=0\text{s}$

$c\left(x,y,z,t_0\right)=f\left(x,y,z\right)\forall \left(x,y,z\right)\in \left\{ℝ,ℝ,{ℝ}^{+}\right\};$ (3)

3) For any from the source, the concentration decreases to zero, i.e:

$\forall x,y,z\text{if}t\to T\text{then}c\to 0;$ (4)

4) The pollutant is not absorbed by the soil surface:

$-K_z\frac{\partial c}{\partial z}=0\text{at}z=0;$ (5)

5) We assume pollutant emissions from a point source located at $\left(0,0,H\right)$ ;

where $f\left(x,y,z\right)=q\delta \left(x\right)\delta \left(y\right)\delta \left(z-H\right)$ ; $q$ is the emission rate and $\delta$ the Dirac’s delta function. So, the mathematical model of pollutant dispersion is presented by the following system:

$\Upsilon =\{\begin{array}{lllll}K\frac{\partial }{\partial t}\left(\nabla c\right)+A\Delta c-{\sigma }^2\nabla c+B\nabla \Delta c\hfill & =\hfill & 0\hfill & \text{in}\hfill & \mathbb{Q},\hfill \\ c\left(x,y,z,0\right)\hfill & =\hfill & f\left(x,y,z\right)\hfill & \text{in}\hfill & {\mathbb{Q}}_0,\hfill \\ K_z\frac{\partial c}{\partial z}\hfill & =\hfill & 0\hfill & \text{at}\hfill & z=0,\hfill \\ c\to 0\hfill & \hfill & \forall x,y,z\hfill & \text{if}\hfill & t\to T,\hfill \\ \left(x_0,y_0,z_0\right)\hfill & =\hfill & \left(0,0,z_0\right)\hfill & \text{if}\hfill & t=0,\hfill \end{array}$ (6)

where $K=\left(\begin{array}{ccc}{K}_x& 0& 0\\ 0& K_y& 0\\ 0& 0& K_z\end{array}\right)$ ; $A=\left(\begin{array}{ccc}{A}_x& 0& 0\\ 0& A_y& 0\\ 0& 0& A_z\end{array}\right)$ ; $B=\left(\begin{array}{ccc}{B}_x& 0& 0\\ 0& B_y& 0\\ 0& 0& B_z\end{array}\right)$ ; ${\sigma }^2=\left(\begin{array}{ccc}{\sigma }_u^2& 0& 0\\ 0& {\sigma }_v^2& 0\\ 0& 0& {\sigma }_w^2\end{array}\right)$ and $A_i=K_i{V}_i$ , $B_i={\sigma }_{\alpha }{\wedge }_i{K}_i$ $\forall i\in \left\{x,y,z\right\}$ and $\alpha \in \left\{U,V,W\right\}$ respectively.

The objective of the mathematical model is to ascertain a solution for the concentration of pollutants, with a view to controlling their dispersion and distribution in the $\mathbb{Q}$ domain. Solving this system (6) of equations directly proves to be a significant challenge due to the inherent complexities of this type of mathematical model. In our study, we opted for a variational methodology to verify the existence and uniqueness of a solution to the system. This approach will also allow us to determine the optimality conditions, enabling us to find an optimal solution.

Let $\mathcal{J}$ be a function denoted by: $\mathcal{J}$ defined by:

$\mathcal{J}:\left(c,c_0\right)\mapsto \mathcal{J}\left(c,c_0\right)={\displaystyle {\int }_{\mathbb{Q}=\Omega \times \left(0,T\right)}L\left(c\left(t\right),c_0\right)\text{d}x\text{d}y\text{d}z\text{d}t},$ (7)

$c\in V_{ad}=\left\{c\in L^{\infty }\left(\mathbb{Q}\right):c_a\le c\le c_b\right\},$

where $c_0$ (initial concentration): the control; $c$ (instantaneous concentration): the state variable; $V_{ad}$ : the control domain; $L$ is the Lagrangian function that meets the following criteria:

• $L:\mathbb{Q}\to ℝ$ is caratheodory function;

• $\forall \left(x,y,z,t\right)\in \mathbb{Q}$ , $L\left(c,c_0\right)$ is of class $C^2\in {ℝ}^2$ ;

• $\forall \left(x,y,z,t\right)\in \mathbb{Q}$ , $L\left(c,c_0\right)$ is convex on the set in $V_{ad}$ .

The problem posed is therefore to minimize this functional $\mathcal{J}\left(c,c_0\right)$ under the system (6) constraint when $c_0\in V_{ad}$ . This is a major constraint of the problem posed. A weak formulation of the problem is formulated by minimizing a new functional $\mathcal{A}$ dependent on $\mathcal{J}\left(c,c_0\right)$ .

Let $\mathcal{A}$ :

$\mathcal{A}\left(c,c_0;p\right)=\mathcal{J}\left(c,c_0\right)+{\displaystyle {\int }_{\mathbb{Q}}ℳ\left(c,c_0\right)\cdot p\text{d}w}$ (8)

$\mathcal{A}\left(c,c_0;p\right)={\displaystyle {\int }_{\mathbb{Q}=\Omega \times \left(0,T\right)}L\left(c\left(t\right),c_0\right)\text{d}w}+{\displaystyle {\int }_{\mathbb{Q}=\Omega \times \left(0,T\right)}ℳ\left(c,c_0\right)\cdot p\text{d}w},$ (9)

with $p\in {\left[C_0^2\left(\mathbb{Q}\right)\right]}^3$ the adjointe variable and we denote $\text{d}x\text{d}y\text{d}z\text{d}t=\text{d}w$ . The function $ℳ\left(c,c_0\right)$ represents the function associated with the variational equation of the evolution equation of the mathematical model of pollutant dispersion.

Furthermore, optimized functionality $\mathcal{J}$ returns to optimized $\mathcal{A}$ :

$\underset{c;c_0\in V_{ad}}{\min }\mathcal{J}\left(c,c_0\right)=\underset{c;c_0}{\min }\underset{p}{\max }\mathcal{A}\left(c,c_0,p\right)$ (10)

Existence and uniqueness

In this session, we use monotone operator theory [22] to demonstrate the existence and uniqueness of a solution to the problem at hand. Operator theory is a powerful tool used in optimization.

In the theoretical study of the model, matrices $A$ , $K$ , ${\sigma }^2$ and $B$ are considered constants.

We consider the evolution equation:

$K\frac{\text{d}}{\text{d}t}\left(\nabla c\right)+A\Delta c-{\sigma }^2\nabla c+B\nabla \Delta c=0\text{in}\mathbb{Q}.$ (11)

Let $ℒ$ the operator associated of the equation (11) such that:

$ℒ(.)=K\frac{\text{d}}{\text{d}t}\left[\nabla (.)\right]+A\Delta (.)-{\sigma }^2\nabla (.)+B\nabla \Delta (.)$ (12)

and

$\begin{array}{l}\mathcal{D}\left(ℒ\right)=\{\varsigma \in H_0^1:ℒ\left(\varsigma \right)=K\frac{\text{d}}{\text{d}t}\left[\nabla \left(\varsigma \right)\right]+A\Delta \left(\varsigma \right)-{\sigma }^2\nabla \left(\varsigma \right)\\ +B\nabla \Delta \left(\varsigma \right);\varsigma \left(x,y,z,0\right)={\varsigma }_0{K}_z\frac{\partial \varsigma }{\partial z}=0\text{if}z=0;\\ \stackrel{}{}\text{and}\varsigma \left(x,y,z,\infty \right)=0\forall \left(x,y,z\right)\in \Omega \}.\end{array}$(13)

Let $c_1$ and $c_2$ , two molarities of the same pollutant. Let’s calculate the dot product in Sobolev space [20]:

$\begin{array}{l}\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \\ =K{\displaystyle {\int }_{\mathbb{Q}}\frac{\text{d}}{\text{d}t}\left[\nabla \left(c_1-c_2\right)\right]\left(c_1-c_2\right)\text{d}w}+A{\displaystyle {\int }_{\mathbb{Q}}\left[\Delta \left(c_1-c_2\right)\right]\left(c_1-c_2\right)\text{d}w}\\ -{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \left(c_1-c_2\right)\right]\left(c_1-c_2\right)\text{d}w}+B{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \Delta \left(c_1-c_2\right)\right]\left(c_1-c_2\right)\text{d}w}.\end{array}$ (14)

Once the calculation and minimization (proof of the calculation and minimization can be found in Appendix A) of this quantity $\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle$ have been carried out, we obtain the following relationship:

$\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \ge {\lambda }^{\prime }{\Vert c_1-c_2\Vert }^2,$ (15)

with ${\lambda }^{\prime }=constant$ . So the operator $ℒ$ is strongly monotone.

Let $c_0$ and $c$ be two different concentrations:

$\begin{array}{c}ℒ\left(c\right)-ℒ\left(c_0\right)=K{\displaystyle {\int }_{\mathbb{Q}}\frac{\text{d}}{\text{d}t}\left[\nabla \left(c-c_0\right)\right]\text{d}w}+A{\displaystyle {\int }_{\mathbb{Q}}\left[\Delta \left(c-c_0\right)\right]\text{d}w}\\ -{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \left(c_1-c_2\right)\right]\text{d}w}+B{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \Delta \left(c-c_0\right)\right]\text{d}w}.\end{array}$ (16)

If $c\to c_0$ therefore $\left(ℒ\left(c\right)-ℒ\left(c_0\right)\right)\to 0$ .

This operator is continue. Thus, the problem (10) has a unique solution.

Optimality systems

We determine the variation of the Lagrange functional and apply the principle of least actions to it in order to find the optimality condition.

Let’s $\delta c$ and $\delta c_0$ be small variations of $c$ and $c_0$ ($c,c_0\in H\left(Q\right)$ ); with $L_c=\frac{\partial L}{\partial c}$ and $L_{c_0}=\frac{\partial L}{\partial c_0}$ .

The variation of Lagrange functional is:

$\begin{array}{c}\Delta \mathcal{A}={\displaystyle {\int }_{\mathbb{Q}}\left[L_c\left(c,c_0\right)-K\nabla \frac{\partial p}{\partial t}+A\Delta p+{\sigma }^2\nabla p+B\Delta \left(\nabla p\right)\right]\delta c\text{d}w}\\ +o\left(\sqrt{{\left(\delta c\right)}^2+{\left(\delta c_0\right)}^2}\right)+{\displaystyle {\int }_{\Omega }\left[{\displaystyle {\int }_0^T{L}_{c_0}\left(c,c_0\right)\text{d}t}+K\nabla p\left(0\right)\right]\delta c_0\text{d}\Omega }\end{array}$ (17)

Proof: See Appendix B for details.

Based on the principles of least action and setting the derivative of the Lagrangian with respect to the variables $c$ and $c_0$ equal to 0, we obtain the following optimality system:

• direct model

$\{\begin{array}{lllll}K\frac{\partial }{\partial t}\left(\nabla c\right)+A\Delta c-{\sigma }^2\nabla c+B\nabla \Delta c\hfill & =\hfill & 0\hfill & \text{in}\hfill & \mathbb{Q},\hfill \\ c\left(x,y,z,0\right)\hfill & =\hfill & f\left(x,y,z\right)\hfill & \text{in}\hfill & {\mathbb{Q}}_0,\hfill \\ K_z\frac{\partial c}{\partial z}\hfill & =\hfill & 0\hfill & \text{at}\hfill & z=0,\hfill \\ c\to 0\hfill & \hfill & \forall x,y,z\hfill & \text{if}\hfill & t\to \pm \infty ,\hfill \\ \left(x_0,y_0,z_0\right)\hfill & =\hfill & \left(0,0,H\right)\hfill & \text{if}\hfill & t=0,\hfill \end{array}$ (18)

• adjoint model

$K\nabla \frac{\partial p}{\partial t}-A\Delta p-{\sigma }^2\nabla p-B\Delta \left(\nabla p\right)=0,p\left(T\right)=0,$ (19)

• optimality condition

$K\nabla p\left(0\right)=0.$ (20)

Analysis of the resulting optimality system leads to the conclusion that an optimal solution exists that minimizes the functional (10). In this study, an analytical solution to the mathematical model will be obtained through the Fourier transform. The analytical (see Appendix C for details) approach is:

$\begin{array}{c}c\left(x,y,z,t\right)=\frac{q}{8\overline{u}{\left(\pi t\right)}^{3/2}}\times \sqrt{\frac{K_x{K}_y{K}_z}{B_x{B}_y{B}_z}}\times \exp [\frac{{\left(x{K}_x+A_{x}t+{\sigma }_u^{2}t\right)}^2}{\left(4B_{x}t\right)}\\ +\frac{{\left(y{K}_y+A_{y}t+{\sigma }_v^{2}t\right)}^2}{\left(4B_{y}t\right)}+\frac{{\left(K_z\left(z\right)+A_{z}t+{\sigma }_w^{2}t\right)}^2}{4B_{z}t}].\end{array}$(21)

2.3.2. Validation of the Model

Parametrization of diffusion coefficients

The implementation of this analytical solution (21) requires the parameterization of diffusion coefficients and standard deviations of velocities. This parameterization depends on the type of stability exhibited by the atmosphere in which the dynamics of the pollutant are studied. There are three types of atmospheric stability: a neutral atmosphere, a stable atmosphere, and an unstable atmosphere. This stability depends on the velocity, the vertical temperature gradient, which in turn depends on the sunshine (day) and the cloud cover (night). Stability is a complex phenomenon of turbulence. Turbulence is represented by a class that allows parameterization of the diffusion standard deviations [28].

BATCHELOR (1949) [29], proposes a parameterization of the diffusion coefficients according to the standard deviations of the concentrations, which results in the following relationships:

$K_{\alpha }=\frac{\overline{u}}{2x}{\sigma }_{\alpha }^2,$ (22)

with $\alpha =\left\{x,y,z\right\}$ , where $\overline{u}$ mean-wind speed; $x$ is downwind distance and ${\sigma }_{\alpha }$ is standard deviation or dispersion parameter in the $\alpha th$ direction. For stable, unstable or neutral conditions, the expression for ${\sigma }_x$ and ${\sigma }_y$ are (CIRILLIO and POLI 1992) [30]:

${\sigma }_x=x{\left[\cosh \left({\sigma }_{\theta }^2\right)-1\right]}^{-1/2},$ (23)

${\sigma }_y=x{\left[\sinh \left({\sigma }_{\theta }^2\right)\right]}^{-1/2},$ (24)

where ${\sigma }_{\theta }$ is the horizontal direction of the wind. The variable ${\sigma }_z$ is parameterized according to the Pasquill-Turner classification table [31]. The diffusion coefficients, $K_x$ and $K_y$ , vary depending on wind speed and horizontal direction, as well as space. Wind speed is a time-dependent variable. Therefore, these coefficients vary over time. The diffusion coefficient $K_z$ varies depending on atmospheric stability and space.

The standard deviations of the velocities are expressed using the following formulas:

• for unstable condition:

${\left(\frac{{\sigma }_u}{\overline{u}}\right)}^2={\left(\frac{{\sigma }_v}{\overline{u}}\right)}^2=0.36{\left(w_{\text{*}}/\overline{u}\right)}^2,$(25)

${\left(\frac{{\sigma }_w}{\overline{u}}\right)}^2=0.16{\left(w_{\text{*}}/\overline{u}\right)}^2,$(26)

• for stable condition:

${\left(\frac{{\sigma }_u}{\overline{u}}\right)}^2=6.25{\left(u_{\text{*}}/\overline{u}\right)}^2,$(27)

${\left(\frac{{\sigma }_v}{\overline{u}}\right)}^2=3.61{\left(u_{\text{*}}/\overline{u}\right)}^2,$(28)

${\left(\frac{{\sigma }_w}{\overline{u}}\right)}^2=1.69{\left(u_{\text{*}}/\overline{u}\right)}^2,$(29)

where $u_{\text{*}}$ is friction velocity and $w_{\text{*}}$ convective velocity. Turbulence intensities depend on average wind speed, convection, friction velocity, and atmospheric stability. Average wind speed, convection, and friction velocity vary over time. Therefore, these intensities vary over time. The necessary and sufficient condition of these empirical relationships is the hypothesis of weak winds in the lower troposphere (atmospheric boundary layer) [32]. Weak winds are observed almost everywhere precisely in the atmospheric boundary layer. Tropical areas are the most exposed to weak winds which promote the diffusion of pollutants [33].

Assimilation data (4D-Var)

A model, however efficient, is not sufficient on its own to make a prediction: it is essential to provide an initial condition deduced from the dynamic observation of the flow [34]. It is essential that this initial condition be as close as possible to the observations while remaining compatible with the model. There are two methods for determining this initial condition: optimal interpolation (OI), which is a statistical method, and the variational method. Optimal interpolation encountered difficulties in the 1980s, notably due to its inability to dynamically handle uncertainties in forecasting, its difficulty in determining initial conditions, and its failure to adequately integrate satellite data [35]. The variational approach 4D-Var (from control theory [36]) formulates data assimilation as a nonlinear optimization problem constrained by a numerical model [37]. Its objective is to take into account recent atmospheric observations in order to readjust the trajectory of a numerical weather prediction model [38]. This approach involves minimizing a functional using the adjoint method in order to find its gradient (Lions 1968 [39]). It found its application in meteorology for the first time by Le Dimet and Talagrand (1986) [40], to one-dimensional models, to filtered models (Lewis et al. 1985) [41]; (Derber 1987) [42], and shallow water equations (Courtier and Talagrand 1990) [43].

However, we use this technique to calibrate our model and determine initial concentrations, which will serve as reference data for forecasts. Thus, the pollutant concentration is evaluated according to time, pressure level, and within each subgrid (cell) of the domain $\mathbb{Q}$ , using the analytical formula (21).

Let $C_0=\left\{c_{0_1},\cdots ,c_{0_N}\right\}$ be the set of initial concentrations evaluated at an altitude and $C_k$ the one evaluated at $t_k$ . Let:

$C_k={ℐ}_{t_0\to t_k}\left(C_0\right)\text{with}{C}_k\in {ℝ}^{d_x};k=1,\cdots ,N$ (30)

where $ℐ:{ℝ}^{d_x}\times {ℝ}^{d_x}$ represents the model describing the evolution of the concentration from $t_0$ to $t_k$ .

$d_x=n_t\times n_{level}\times n_{lat}\times n_{lon}$ (31)

$d_x=8\times 5\times 88\times 121$

$d_x=425920,$ (32)

$n_t=8$ is the number of hour steps in a day; $n_{level}=5$ is the number of pressure levels; $n_{lat}=88$ and $n_{lon}=121$ are the dimensions of the domain where the cells of dimensions are 28 km × 28 km. Let:

$\mathscr{H}:{ℝ}^{d_x}\to {ℝ}^{d_y},$ (33)

be the observation operator that links the estimated state of the model to the observation corresponding to each instant $t_k$ with $k=1,\cdots ,N$ and $d_y=5\times 8\times 30\times 41=49200$ . In mathematical modeling, the observation operator $\mathscr{H}$ maps the complete state vector of the four-dimensional model (time × pressure levels × latitude × longitude) into the observation space. This representation is achieved through a linear projection. The model state is first vectorized into a one-dimensional array: $C_k\in {ℝ}^{d_x}$ . The operator $\mathscr{H}$ is structured as a sparse linear selection matrix, whose function is to extract only the state components corresponding to the available satellite observations. More specifically, a Boolean observation mask is used to identify the observed grid points and pressure levels. The operator, by its linearity, is a function of time, as demonstrated by the observation mask, which reflects the spatial and temporal availability of satellite retrievals.

The observation vector $y_k^0$ at time $t_k$ and the noise vector ${\nu }_k$ are defined by the following application:

$y_k^0=\mathscr{H}\left(C_k\right)+{\nu }_k,$ (34)

with $y_k^0\in {ℝ}^{d_y}$ ; ${\nu }_k~\mathcal{N}\left({\sigma }^2,{ℛ}_k\right)$ and ${\sigma }^2={0.01}^2$ [24] [37]. The observation error (${\nu }_k$ ) is unbiased, uncorrelated in time. The observation vector ($y_k^0$ ) corresponds to satellite images of carbon monoxide (CO) concentration at a time $t_k$ . The functional associated with the variational approach is [44]:

$\begin{array}{c}\mathbb{S}\left(C_0\right)=\frac{1}{2}{\displaystyle \sum }_{k=0}^N{\left[y_k^0-\mathscr{H}\left(C_k\right)\right]}^{\text{T}}{R}_k^{-1}\left[y_k^0-\mathscr{H}\left(C_k\right)\right]\\ +\frac{1}{2}{\left[C_0-C_b\right]}^{\text{T}}{ℱ}_0^{-1}\left[C_0-C_b\right],\end{array}$ (35)

with

• $C_b~\mathcal{N}\left(0.04,{ℱ}_0\right)$ is the initial background state [24];

• $R_k={0.01}^2{I}_{dy}\in {ℝ}^{d_y\times d_y}$ is the error covariance diagonal matrix of the observation [24] [37];

• ${ℱ}_0=0.04I_{d_x}\in {ℝ}^{d_x\times d_x}$ the model error covariance diagonal matrix [37],

with $I_{d_x}$ the identity matrix.

The problem of data assimilation using the 4D-Var method is therefore:

$\left(ℐ\mathcal{V}\right):\{\begin{array}{ll}{C}_k={ℐ}_{t_0\to t_k}\left({\widehat{C}}_0\right)\hfill & C_k\in {ℝ}^{d_x},\hfill \\ y_k^0=\mathscr{H}\left(C_k\right)+{\nu }_k\hfill & y_k^0\in {ℝ}^{d_y},\hfill \\ C\left(0\right)=C_0\hfill & C_0\in {ℝ}^{d_x},\hfill \\ {\widehat{C}}_0=\arg \underset{C_0\in {ℝ}^{d_x}}{\min }\mathcal{S}\left(C_0\right)\hfill & {\widehat{C}}_0\in {ℝ}^{d_x}.\hfill \end{array}$ (36)

The control vector consists of the complete four-dimensional concentration field defined over all time steps of the assimilation window, all pressure levels, and all horizontal grid points. The algorithm of choice for minimisation is the Quasi-Newton limited memory L-BFGS algorithm as devised in (Nocedal 1980) [45] and (Liu and Nocedal 1989) [46]. The 4D-Var assimilation window spans a period of three hours (03 h), corresponding to the complete integration period of the direct model. The model is integrated continuously over this fixed window. Consequently, all satellite observations available during this period are assimilated simultaneously at their respective acquisition times. The 4D-Var cost function is minimized by the quasi-Newton L-BFGS algorithm with an analytically derived adjoint gradient. It should be noted that the number of iterations is limited to 100 ($maxiter=100$ ), and that the default stopping criteria of the SciPy implementation are applied (projected gradient norm and relative cost reduction). The corresponding inverse covariance matrices weight the observation and background terms of the cost function and ensure regularization of the problem. The dynamic system consists of nine (09) input data points. These include eight (8) ERA5 meteorological variables, as well as carbon monoxide concentration images. The data are daily with a three-hourly interval. The model makes eight (08) estimates per day based on altitude. Due to the considerable size of the data, the numerical simulations require the use of a powerful and suitable computer, given the complexity of the necessary mathematical calculations. Due to the size of the available data, the simulations were performed on a calculator. To evaluate the robustness of the model, we use the Root Mean Square Error (RMSE) metric [24], the fractional bias (FB) [23], and the fraction of predictions within a factor of two of the observations (FAC2) [23]:

$RMSE=\sqrt{\frac{1}{M}\cdot {\displaystyle \sum }_{k=0}^M{\left({\widehat{c}}_{0_k}-y_k^0\right)}^{\text{T}}\cdot \left({\widehat{c}}_{0_k}-y_k^0\right)},$(37)

$FB=\frac{\overline{y_k^0}-\overline{{\widehat{c}}_{0_k}}}{0.5\left(\overline{y_k^0}+\overline{{\widehat{c}}_{0_k}}\right)},$ (38)

$FAC2=\frac{{\widehat{c}}_{0_k}}{y_k^0},$ (39)

where $\overline{C}$ represents the average over the dataset and ${\widehat{c}}_{0_k}$ the estimated initial concentration at time $t_k$ .

3. Results and Discussion

In this session, we present the results of the implementation of the model coupled with the variational approach (4D-Var), which is divided into two parts. First, we compare the output maps of the coupled model with the observations incorporated into it. The second part of this study focuses on the analysis of the output maps of the coupled system. It should be noted that the experiments were conducted using data from the dry season (December-January-February), specifically December 1^st, 2023. The dry season is characterized by lower humidity and higher than normal temperatures along the West African coast. The prevailing wind during this season carries many particles from the Sahara Desert southward. This is a period during which high levels of carbon monoxide (CO) accumulate in the Gulf of Guinea [13] [47]. The weather conditions on the chosen day are determined by solar radiation during the day and cloud cover at night, according to the Pasquill-Turner classification table [31].

3.1. Algorithm Workflow: Coupled System Algorithm Model-4D-Var

Figure 2 shows satellite images (left) and dispersion maps (right) from the Model-4D-Var coupling system. These dispersion maps illustrate the evolution of dispersion and the estimated concentration of carbon monoxide at 1000 hPa, 950 hPa, and 850 hPa as a function of time. We observe that the dispersion maps of the coupling system (right-hand images) are smoother and allow for the quantification of CO pollutant plumes. Furthermore, these maps enable us to track the spatio-temporal trajectories of the plumes with greater precision. They contribute to a better understanding of pollutant transport at different meteorological scales in the Gulf of Guinea. This is due to the integration of the 4D-Var data assimilation technique, a major strength of the system. This technique allows it to produce better estimates consistent with satellite observations and to reduce uncertainties compared to available emission inventories, which are sometimes biased in West Africa. This performance is close to modern techniques based on the automation of quantification and detection of plumes of pollutants (CH₄) recently developed [48]-[52]. These machine learning techniques use satellite and hyperspectral images (Sentinel-2, PRISMA, EnMAP, EMIT and CarbonMapper) to detect plumes and estimate concentrations of pollutants from combustion and industries etc.

Obtaining carbon monoxide dispersion maps with a resolution scale of (0.25˚ × 0.25˚) compared to satellite images with a resolution of 0.75˚ × 0.75˚ is part of one of the model’s performances. The concentration in each grid cell is estimated using the analytical formula (21). This formula takes into account the horizontal

Figure 2. Satellite images are on the left and those generated by the model are on the right at 1000 hPa, 950 hPa 850 hPa depending on the time.

and vertical distribution of carbon monoxide as a function of pressure levels. It considers advection, diffusion, convection, atmospheric stability (clouds, solar radiation), and turbulence through parameterizations of the diffusion coefficients $\left(k_x,K_y,K_z\right)$ and velocity standard deviations $\left({\sigma }_u,{\sigma }_v,{\sigma }_w\right)$ .

The low values of the mean squared error (Figure 3), attest that the results of the model are close to the observations. This property also demonstrates the system’s effectiveness and its ability to generate dispersion maps that accurately reflect real-world observations. The system reproduces carbon monoxide dispersion more accurately than observations. Indeed, at 1000 hPa these values are considered high compared to other pressure levels, namely (950 hPa; 900 hPa; 850; 800 hPa) between 00 h - 06 h. The results obtained are inconsistent with observations recorded at this pressure level 1000 hPa. The anomaly observed at 1000 hPa is related to the choice of parameters for the diffusion coefficients and velocity standard deviations. At this level, the appropriate parameters for the diffusion coefficients are those of Monin and Obukhov [53]. Those suitable for the velocity standard deviations are those proposed by Irwin (1970) and Wyngaard et al. (1974-1975), as reported by Hanna et al. [27].

Figure 3. RMSE, fractional bias and FAC2 of the model.

However, the objective of our study is to investigate the dispersion of pollutants at pressure levels above 1000 hPa; therefore, we cannot use these formulations, which could introduce physical inconsistencies into the model. It should be noted, however, that this anomaly does not significantly affect the overall performance of the model when the pressure level is above 1000 hPa. The fractional bias (FB) values obtained at each pressure level were found to be very low and negative ($FB<0$ ). The negative sign ($FB<0\Rightarrow \overline{y_k^0}<\overline{c_{0_k}}$ ) indicates that the model overestimates carbon monoxide concentration values. However, the results obtained indicate that the overestimates is negligible compared to the actual values from satellite images. This study shows that the model estimates correlate perfectly with satellite images.

When analyzing the data (Table A1), it is essential to specify that the fraction data (FAC2) values, determined for each pressure level (Figure 3), are between 0.5 and 2. These values comply with current standards and reflect the model’s good performance ($0.5\le FAC2\le 2$ ) in terms of reproducing the pollutant dispersion process. In addition, it has been found that the model’s performance is limited to 1000 hPa between midnight and 6 a.m. This observation confirms the anomaly identified by the RMSE at this constant pressure level.

3.2. Analysis of the Dispersion Output Maps of the Model-4D-Var

The 1000 hPa dispersion maps (Figure 4) show carbon monoxide concentration

Figure 4. Temporal evolution of carbon monoxide concentration at 1000 hPa.

values ranging from 2.5 ppmv to 1.4 ppmv (in practice 1 ppmv = 1 ppm), with a resolution of three hours. According to the World Health Organization (WHO), the recommended atmospheric concentration threshold is 7 mg/m³ (approximately 6.1 ppm at 25˚C) for a 24 h exposure period [54]. Further research has established a correlation between exposure to ambient carbon dioxide concentrations below 1 ppm (approximately 1.15 mg/m³ at 25˚C) over a period of one hour and an increased risk of cardiovascular disease [55] [56]. This observation contradicts the guidelines established by the World Health Organization (WHO). Chen et al. (2021) [57] conducted studies using data from 1979 to 2016, covering 337 cities in 18 countries. This research highlights the absence of a threshold concentration for CO, suggesting instead a correlation between the risk of cardiovascular disease or daily mortality and a 1 mg/m³ (approximately 0.87 ppm at 25˚C) increase in this gas. In their study, Samoli et al. (2007) [58] demonstrated the existence of significant co-mortality risks for concentrations below 0.5 mg/m³ (approximately 0.44 ppm at 25˚C), in the case of maximum average exposure over an 8 h period.

These results show that, for exposure lasting one hour (1 h) and a threshold value of 1 ppmv, the risk of cardiovascular disease is proven. Furthermore, for a period of three hours, the value must be below 1 ppmv. Polluted areas are therefore territories in which carbon monoxide concentrations exceed the threshold of 1 ppmv. In these areas, populations are exposed to levels of carbon monoxide (CO) that can have harmful effects on health, including respiratory diseases and even death. It has been found that most of these areas are those where continuous emissions of pollutants occur. These emissions are characterized by their spread over altitude and increased persistence on a satellite scale.

3.2.1. Quantification of Local Plumes

• At 1000 hPa (Figure 5) and at 06 h, we detect a local plume that appears in Liberia in the region of 5.5˚N and 10˚W with a concentration of approximately to 2.25 ppmv. This plume disperses slowly under the influence of a light southerly wind, and disappears around 12 h. It has little impact on the area during this period.

• At 950 hPa (Figure 6), a second plume is located south of Nigeria at around 7˚E and 6˚N. It is present at 00 h with a maximum concentration of approximately 1.6 ppmv, which decreases over time (00 h; 03 h; 06 h), but rises again around 18 h with a value close to 0.45 ppmv.

The dispersal of this plume is primarily influenced by the monsoon (south wind) due to its location near the Atlantic Ocean. This south wind, observed in the southern coastal regions, blows between 00 h and 06 h, and ceasing when the harmattan becomes dominant. Before 06 h, we note a sunrise with a less cloudy sky (Figure 7) marking the dispersion of the plume. During the sunny day, the sky is clear and the daytime heating becomes significant. This marks the presence of the sea breeze which promotes a rapid dispersion of the plume with a decrease in concentration, or even cancellation between 12 h - 03 h (Figure 6). Clouds reappear after 06 h, the atmosphere becomes stable, and dispersion is limited, leading to an accumulation of pollutants in the area. This accumulation pollutes the air and poses a risk to the local population. The presence of carbon monoxide in the southern region of Nigeria has been highlighted by the work of some authors.

Figure 5. Temporal evolution of CO dispersion with interpolation of wind at 1000 hPa (Images generated by the Model-4D-Var).

Figure 6. Temporal evolution of carbon monoxide dispersion with interpolation of wind at 950 hPa (Images generated by the Model-4D-Var).

Figure 7. Illustration of cloud cover fraction and solar radiation.

Gijs et al. (2023) [2]; Mayowa et al. (2024) [3]; Ajoke et al. (2025) [4] highlighted the presence of carbon monoxide in southern and western Nigeria, caused by biomass burning and gas flaring, as well as by forest fires, incomplete combustion in road transport, residential heating and industry, and by transport and urbanization. An emission source located at Warri (5.33˚N, 5.36˚E) is specified by Ajoke et al. (2025) [4] in their work. This source detected at an altitude of 100 m near the ground helps to explain the origin of the plume identified by our model at an altitude of approximately 500m in this area.

3.2.2. Quantification of Continuous Vertical Plumes

Three plumes persisted at all altitudes considered in our study. They exhibited a spatio-temporal distribution with variability in concentration at each altitude.

• A plume identified north of Nigeria, near the city of Kano, whose origin is located near 9.30˚E and 12.30˚N. At 1000 hPa (Figure 8), it is less persistent between 00 h and 03 h but more persistent from 06 h onwards. We note a significant spatial occupation of the plume with a considerable concentration. The plume’s distribution starts in northern Nigeria, crosses northern Benin from 12 h and extends to northern Togo. Starting at 950 hPa, the plume concentrates at its origin at 00 h with a molarity close to 0.6 ppmv. After 00 h, under the effect of a high-intensity dry wind, it disperses rapidly with a decrease in concentration. This means that a significant concentration of carbon monoxide is being emitted at a low altitude in northern Nigeria, and is visible from space. The results of the work carried out on the quantification of carbon monoxide from satellite images of TROPOMI (Sentinel-5 Precusor) by Gijs et al. (2023) [2], illustrate the presence of CO emission sources in the city of Kano (Nigeria) with an emission rate of 0.53 (0.41 - 0.62) Tg.yr⁻¹. The dispersion of the pollutant from these sources confirms the plumes represented by our model.

• A second plume has been detected in Ghana. Analysis of the wind circulation direction made it possible to locate an emission source which is located in the city of Kumasi (6.36˚N; 1.37˚W). This emission is due to transport, industrial

Figure 8. Temporal evolution of carbon monoxide dispersion with interpolation of wind at 1000 hPa and 950 hPa (Images generated by the Model-4D-Var).

• Activities and combustion. Ajoke et al. (2025) [4] indicated in their work, by approximation, the existence of a carbon monoxide emission source close to this city. The results of Gijs et al. (2023) [2] confirm the existence of a carbon monoxide emission source in the city of Kumasi. This result provides further details and confirms the presence of our plumes in this area. The dispersion of the pollutant emitted in this area is influenced by the wind from the Southwest at 1000 hPa. At low altitude (1000 hPa, Figure 5) between 03 h - 09 h, the pollutant disperses from the source, under the effect of the land breeze (south-west wind). It spreads south and extends south into Togo and Benin with a low concentration that varies between 1 ppmv - 1.5 ppmv.

• A third plume which is located in the Central African Republic, at 7.55˚N and 20˚E (Figure 8). It is located in an area characterized by the Sudanese savanna and persists at all altitudes considered in this study. The dispersion of this plume influences the entire Gulf of Guinea area. It comes mainly from forest and agricultural fires [13]. At low altitudes, it is less present at night with a low concentration. However, during the day, a significant concentration of carbon monoxide is observed, which disperses into other countries and contributes to increasing the concentration of other plumes originating from other sources located in its direction. The distribution of this plume is strongly favored by the northeast wind. Onojeghuo et al. (2025) [4] indicated the presence of a source emitting this pollutant in this area.

4. Conclusions and Suggestions

This article involves the development of a theoretical and numerical model to explain the atmospheric dispersion of carbon monoxide in the Gulf of Guinea. Theoretically, the model is based on the scalar flux transport equation. The existence, uniqueness, and optimality conditions of this model are studied using the tools of control theory.

In practice, an analytical expression is derived from this theoretical model. This expression combines several parameters likely to influence the phenomenon of pollutant dispersion in the Gulf of Guinea. This analytical approach is coupled with a 4D-Var filter. The assimilation of ERA5 reanalysis data (horizontal and vertical wind, temperature, cloud cover, convection velocity, solar radiation, friction velocity) and carbon monoxide data from a coarse grid allowed the model to reproduce dispersion maps consistent with satellite images. The output maps have a higher resolution than satellite images. The consideration of atmospheric stratification and the variability of wind types in the Gulf of Guinea demonstrates the model’s ability to adapt to the study area. The results obtained represent a significant contribution to the understanding of carbon monoxide transport in the Gulf of Guinea. Looking ahead, integrating the Kalman filter into the model will improve its future predictions.

Acknowledgements

The authors express their gratitude to the IMSP (Institut de Mathématiques et de Sciences Physiques) for providing its high-performance computing resources, which enabled the numerical simulations to be performed.

This research did not receive any specific funding from public, commercial, or non-profit funding organizations.

Appendix A: Calculation and Minimization of Quantity: $ℒ\left(c_1\right)-ℒ\left(c_2\right)$

Let’s calculate $ℒ\left(c_1\right)-ℒ\left(c_2\right)$ :

$\begin{array}{c}ℒ\left(c_1\right)-ℒ\left(c_2\right)=K\frac{\text{d}}{\text{d}t}\left[\nabla \left(c_1-c_2\right)\right]+A\Delta \left(c_1-c_2\right)-{\sigma }^2\nabla \left(c_1-c_2\right)\\ +B\nabla \Delta \left(c_1-c_2\right).\end{array}$ (40)

Let’s calculate the dot product of this difference (40) in Sobolev space:

$\begin{array}{l}\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \\ =K{\displaystyle {\int }_{\mathbb{Q}}\frac{\text{d}}{\text{d}t}\left[\nabla \left(c_1-c_2\right)\right]\left(c_1-c_2\right)\text{d}w}+A{\displaystyle {\int }_{\mathbb{Q}}\left[\Delta \left(c_1-c_2\right)\right]\left(c_1-c_2\right)\text{d}w}\\ -{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \left(c_1-c_2\right)\right]\left(c_1-c_2\right)\text{d}w}+B{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \Delta \left(c_1-c_2\right)\right]\left(c_1-c_2\right)\text{d}w}.\end{array}$ (41)

Let $\left(c_1-c_2\right)=h\left(t\right)g\left(x,y,z\right)$ , with $h\left(t\right)\in \left[0,T\right]$ and $g\left(x,y,z\right)\in \Omega$ , hence we obtain:

$\begin{array}{l}\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \\ =K{\displaystyle {\int }_{\mathbb{Q}}\frac{\text{d}{h}^2}{\text{d}t}\left[\nabla g\right]g\text{d}w}+A{\displaystyle {\int }_{\mathbb{Q}}\left[\Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\\ -{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \left(hg\right)\right]\left(hg\right)\text{d}w}+B{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\\ =K{\displaystyle {\int }_{\Omega }\left[\nabla g\right]g\text{d}\Omega }{\displaystyle {\int }_0^T\frac{\text{d}{h}^2}{\text{d}t}\text{d}t}+A{\displaystyle {\int }_{\mathbb{Q}}\left[\Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\\ -{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \left(hg\right)\right]\left(hg\right)\text{d}w}+B{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\\ =\frac{K}{2}{\displaystyle {\int }_{\Omega }\left[\nabla g\right]g{\left[h^2\right]}_0^T\text{d}\Omega }+A{\displaystyle {\int }_{\mathbb{Q}}\left[\Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\\ -{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \left(hg\right)\right]\left(hg\right)\text{d}w}+B{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\end{array}$

$\begin{array}{l}\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \\ =\frac{K}{2}{\displaystyle {\int }_{\Omega }{h}^{2}T\left[\nabla g\right]g\text{d}\Omega }+A{\displaystyle {\int }_{\mathbb{Q}}\left[\Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\\ -{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \left(hg\right)\right]\left(hg\right)\text{d}w}+B{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \Delta \left(hg\right)\right]\left(hg\right)\text{d}w}.\end{array}$ (42)

We minimize each quantity of the scalar product (42). We suppose that ${\lambda }_i\in {ℝ}_{+}^{\text{*}}$ , with $i=1,2,\cdots ,n$ :

${\lambda }_1\le \Delta \left(hg\right)\le {\lambda }_2$

${\lambda }_1\left(hg\right)\le \left[\Delta \left(hg\right)\right]\left(hg\right)\le {\lambda }_2\left(hg\right)$

$A{\displaystyle {\int }_{\mathbb{Q}}{\lambda }_1\left(hg\right)\text{d}w}\le A{\displaystyle {\int }_{\mathbb{Q}}\left[\Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\le A{\displaystyle {\int }_{\mathbb{Q}}{\lambda }_2\left(hg\right)\text{d}w}$ (43)

${\lambda }_3\le \left[\nabla \left(hg\right)\right]\le {\lambda }_4$

$\left(hg\right){\lambda }_3\le \left[\nabla \left(hg\right)\right]\left(hg\right)\le {\lambda }_4\left(hg\right)$

${\displaystyle {\int }_{\mathbb{Q}}{\lambda }_3\left(hg\right)\text{d}w}\le {\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \left(hg\right)\right]\left(hg\right)\text{d}w}\le {\displaystyle {\int }_{\mathbb{Q}}{\lambda }_4\left(hg\right)\text{d}w}$

$-{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}{\lambda }_4\left(hg\right)\text{d}w}\le -{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \left(hg\right)\right]\left(hg\right)\text{d}w}\le -{\sigma }^2{\displaystyle {\int }_{\mathbb{Q}}{\lambda }_3\left(hg\right)\text{d}w}$ (44)

${\lambda }_5\le \left[\nabla \text{Δ}\left(hg\right)\right]\le {\lambda }_6$

${\displaystyle {\int }_{\mathbb{Q}}{\lambda }_5\left(hg\right)\text{d}w}\le {\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\le {\displaystyle {\int }_{\mathbb{Q}}{\lambda }_6\left(hg\right)\text{d}w}$

$B{\displaystyle {\int }_{\mathbb{Q}}{\lambda }_5\left(hg\right)\text{d}w}\le B{\displaystyle {\int }_{\mathbb{Q}}\left[\nabla \Delta \left(hg\right)\right]\left(hg\right)\text{d}w}\le B{\displaystyle {\int }_{\mathbb{Q}}{\lambda }_6\left(hg\right)\text{d}w}.$ (45)

The combination of equations (43), (44) and (45) makes it possible to obtain a framework for the scalar product:

$\begin{array}{l}{\displaystyle {\int }_{\mathbb{Q}}\left\{A{\lambda }_1\left(hg\right)-{\sigma }^2{\lambda }_4\left(hg\right)+B{\lambda }_5\left(hg\right)\right\}\text{d}w}+\frac{K}{2}{\displaystyle {\int }_{\Omega }{h}^2\left(T\right)\left[\nabla g\right]g\text{d}\Omega }\\ \le \langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \\ \le {\displaystyle {\int }_{\mathbb{Q}}\left\{A{\lambda }_2\left(hg\right)-{\sigma }^2{\lambda }_3\left(hg\right)+B{\lambda }_6\left(hg\right)\right\}\text{d}w}+\frac{K}{2}{\displaystyle {\int }_{\Omega }{h}^2\left(T\right)\left[\nabla g\right]g\text{d}\Omega },\end{array}$ (46)

hence

$\begin{array}{l}\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \\ \ge {\displaystyle {\int }_{\mathbb{Q}}\left\{A{\lambda }_1\left(hg\right)-{\sigma }^2{\lambda }_4\left(hg\right)+B{\lambda }_5\left(hg\right)\right\}\text{d}w}+\frac{K}{2}{\displaystyle {\int }_{\Omega }{h}^2\left(T\right)\left[\nabla g\right]g\text{d}\Omega }\end{array}$ (47)

$\begin{array}{l}\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \\ \ge {\displaystyle {\int }_{\mathbb{Q}}\left\{A{\lambda }_1-{\sigma }^2{\lambda }_4+B{\lambda }_5\right\}\left(hg\right)\text{d}w}+\frac{K}{2}{\displaystyle {\int }_{\Omega }{h}^2\left(T\right)\left[\nabla g\right]g\text{d}\Omega }\end{array}$ (48)

$\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \ge D{\displaystyle {\int }_{\mathbb{Q}}\left(hg\right)\text{d}w}+K{\displaystyle {\int }_{\mathbb{Q}}\frac{\text{d}}{\text{d}t}\left[\nabla \left(hg\right)\right]\left(hg\right)\text{d}w}$ (49)

$\begin{array}{l}\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \\ \ge D{\displaystyle {\int }_{\mathbb{Q}}\left(c_1-c_2\right)\text{d}w}+K{\displaystyle {\int }_{\mathbb{Q}}\frac{\text{d}}{\text{d}t}\left[\nabla \left(c_1-c_2\right)\right]\left(c_1-c_2\right)\text{d}w}\end{array}$ (50)

$\begin{array}{l}\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \\ \ge \min \left(D,K\right){\displaystyle {\int }_{\mathbb{Q}}\left\{\frac{\text{d}}{\text{d}t}\left[\nabla \left(c_1-c_2\right)\right]\left(c_1-c_2\right)+\left(c_1-c_2\right)\right\}\text{d}w}\end{array}$ (51)

$\langle ℒ\left(c_1\right)-ℒ\left(c_2\right),\left(c_1-c_2\right)\rangle \ge {\lambda }_7{\Vert c_1-c_2\Vert }^2,$ (52)

with $D=A{\lambda }_1-{\sigma }^2{\lambda }_4+B{\lambda }_5$ and ${\lambda }_7=\min \left(D,K\right)$ .

Appendix B: Variation of Lagrange Functional

Let’s $p\in {\left[C_0^2\left(Q\right)\right]}^3$ the adjoint variable and $c\in H\left(Q\right)$ .

Remark: ${\left[C_0^2\left(Q\right)\right]}^3=C_0^2\left(Q\right)\times C_0^2\left(Q\right)\times C_0^2\left(Q\right)$ . Let’s $\delta c$ and $\delta c_0$ be small variations of $c$ and $c_0$ .

$\mathcal{A}\left(c,c_0;p\right)={\displaystyle {\int }_{\mathbb{Q}}L\left(c\left(t\right),c_0\right)\text{d}w}+{\displaystyle {\int }_{\mathbb{Q}}ℳ\left(c,c_0\right)\cdot p\text{d}w}$ (53)

$\begin{array}{l}\mathcal{A}\left(c+\delta c,c_0+\delta c_0;p\right)\\ ={\displaystyle {\int }_{\mathbb{Q}}L\left(c+\delta c,c_0+\delta c_0\right)\text{d}w}+{\displaystyle {\int }_{\mathbb{Q}}ℳ\left(c+\delta c,c_0+\delta c_0\right)\cdot p\text{d}w}\end{array}$ (54)

The variation of Lagrange functional is:

$\text{Δ}\mathcal{A}=\mathcal{A}\left(c+\delta c,c_0+\delta c_0;p\right)-\mathcal{A}\left(c,c_0;p\right)$ (55)

$\begin{array}{c}={\displaystyle {\int }_{\mathbb{Q}}L\left(c+\delta c,c_0+\delta c_0\right)\text{d}w}+{\displaystyle {\int }_{\mathbb{Q}}ℳ\left(c+\delta c,c_0+\delta c_0\right)\cdot p\text{d}w}\\ -{\displaystyle {\int }_{\mathbb{Q}}L\left(c,c_0\right)\text{d}w}-{\displaystyle {\int }_{\mathbb{Q}}ℳ\left(c,c_0\right)\cdot p\text{d}w}\end{array}$ (56)

$\begin{array}{c}={\displaystyle {\int }_{\mathbb{Q}}\left\{L\left(c+\delta c,c_0+\delta c_0\right)-L\left(c,c_0\right)\right\}\text{d}w}\\ +{\displaystyle {\int }_{\mathbb{Q}}\left\{ℳ\left(c+\delta c,c_0+\delta c_0\right)-ℳ\left(c,c_0\right)\right\}\cdot p\text{d}w}\end{array}$ (57)

$\Delta \mathcal{A}=\Delta {\mathcal{A}}_1+\Delta {\mathcal{A}}_2,$ (58)

with

$\text{Δ}{\mathcal{A}}_1={\displaystyle {\int }_{\mathbb{Q}}\left\{L\left(c+\delta c,c_0+\delta c_0\right)-L\left(c,c_0\right)\right\}\text{d}w},$ (59)

$\Delta {\mathcal{A}}_2={\displaystyle {\int }_{\mathbb{Q}}\left\{ℳ\left(c+\delta c,c_0+\delta c_0\right)-ℳ\left(c,c_0\right)\right\}\cdot p\text{d}w}.$ (60)

Lemma [20]: Let $c_0\in L^{\infty }\left(\mathbb{Q}\right)$ be the optimal control of the problem and $c$ the corresponding state.

Then, the functional ${\mathcal{A}}_1$ is Frechet differentiable and the variation is expressed as:

$\text{Δ}{\mathcal{A}}_1={\displaystyle {\int }_{\mathbb{Q}}\left[L_c\left(c,c_0\right)\delta c+L_{c_0}\left(c,c_0\right)\delta c_0\right]\text{d}w}+o\left(\sqrt{\delta c^2+\delta c_0^2}\right),$ (61)

where $L_c$ and $L_{c_0}$ stand for the derivative with respect to c and $c_0$ , respectively.

However, let us determine $\text{Δ}{\mathcal{A}}_2$ .

$ℳ\left(c,c_0\right)=K\frac{\partial }{\partial t}\left(\nabla c\right)+A\Delta c-{\sigma }^2\nabla c+B\nabla \Delta c$ (62)

$\begin{array}{l}ℳ\left(c+\delta c,c_0+\delta c_0\right)\\ =K\frac{\partial }{\partial t}\nabla \left(c+\delta c\right)+A\Delta \left(c+\delta c\right)-{\sigma }^2\nabla \left(c+\delta c\right)+B\nabla \Delta \left(c+\delta c\right)\end{array}$ (63)

$\begin{array}{l}ℳ\left(c+\delta c,c_0+\delta c_0\right)-ℳ\left(c,c_0\right)\\ =K\frac{\partial }{\partial t}\left(\nabla \delta c\right)+A\Delta \delta c-{\sigma }^2\nabla \delta c+B\nabla \Delta \delta c.\end{array}$ (64)