Oil Discharge Trajectory Simulation at Selected Baltic Sea Waterway under Variability of Hydro-Meteorological Conditions

Abstract

The paper deals with an important issue related to the identification, modelling, and prediction of environmental pollution in aquatic ecosystems of the Baltic Sea caused by anthropopressure. Water ecosystems are in danger nowadays because of the negative influence of chemical releases in seas, oceans, or inland waters. The crucial issue is to prevent the oil spills and mitigate their consequences. Thus, there is a need for methods capable of reducing the water pollution and enhancing the effectiveness of port and marine environment preservation. The challenge in implementing actions to remove and prevent horizontal oil discharge lies in accurately determining its shape and direction of oil spreading. The author employed a self-designed software utilizing modified and developed mathematical probabilistic models to forecast the movement and dispersion of an oil spill in diverse hydrological and meteorological conditions. This involved determining the trajectory and movement of a spill domain, which consists of elliptical sub-domains undergoing temporal changes. The research results obtained are the initial results in the oil spill simulation problem. This approach represents an expanded and innovative method for determining the spill domain and tracking its movement, applicable to oceans and seas worldwide. It expands upon the methodologies firstly discussed, thereby broadening the range of available techniques in this field. A simple model of an oil spill trajectory simulation and a surface oil slick as an ellipse is illustrated using a time-series of selected hydro-meteorological factors that change at random times. The author proposes a Monte Carlo simulation method to determine the extent of an oil spill in an aquatic ecosystem, taking into account the influence of varying hydro-meteorological conditions. A semi-Markov model is defined to capture the dynamics of these conditions within the spill area and develop an enhanced algorithm for predicting changes in the shape and movement of the spill domain under changing these conditions. By applying the algorithm, a simulation is conducted to provide short-term prediction of the oil discharge trajectory in a selected Baltic Sea waterway. To enhance the accuracy of predicting the process of changing conditions, uniformly tested joint datasets from the open sea water area were incorporated. Finally, the potential future prospects and directions for further research in this field are discussed.

1. Introduction

In recent years, a large threat of chemical releases into international waters has been observed. This pollution means that the concentration of toxic substances is harmful to the environment and local economy [1,2], which implies that the effects can harm biological life and affect human activities. Any intentional or accidental action that results in the spill, emptying, dumping, pumping, emitting, or pouring of oil into water is considered a discharge of oil [3]. The main causes of pollution are seabed exploration and drilling platform construction, oil transport by ships and pipelines, pollution during drilling and pipelines operation, contamination from rivers, and illegal oil discharges from ships [4]. There were four medium spills (7–700 tonnes) and three large spills (>700 tonnes, in Asia and Africa) in 2022 with total volume of around 15,000 tonnes [5]. The current trends seem to be downward, however some fluctuations still occur (Figure 1 and Figure 2).

The significant issue is to forestall the oil slicks and mitigate their cost to the nations and wildlife affected. The oil spill on the water surface behaves very differently depending on its chemical composition as well as on the hydrological parameters of the basin [1]. The removal of pollution, especially oil discharge in the open sea area, is very difficult due to the variability of hydro-meteorological conditions [1,2]. The spreading of oil generally depends on wind, wave height, and currents [1,6].

In terms of transport intensification, the Baltic Sea is widely recognized as one of the most congested seas globally. In the era of globalization of the economy, we can observe a negative impact of chemical releases, especially during winter, when the sun light falls on the Earth with less intensity and the days are shorter [7,8,9,10]. In terms of threats to the Baltic Sea, maritime transport is one of the most dangerous factors affecting the environment. That is why safety issues should be considered, especially in the Baltic Proper basin [1]. Difficulties associated with reaching the scene of the event, carrying out the action, the risk of sinking the ship, including fuel sunk and large amounts of dangerous cargo on board, cause disasters on a scale of much larger areas compared to disasters on land [1]. Due to the amount of transported cargo, the greatest threats of an ecological disaster as a result of accidents come from oil tankers and container ships. It will certainly not be indifferent to the flora and fauna of the sea, and thus to other sea users. There is also the fear that the effects of contamination will be felt for a long time after the accident occurred.

This research problem of sea pollution identification, modelling, and prediction due to the spills of hazardous substances has existed for many years and is becoming more and more significant due to the intensive development of marine water purity monitoring systems [5,6,7,8,9]. The increase in computerization, automation, and miniaturization of measurement and navigation systems leads to the prevention and quick detection of oil spills by means of monitoring activities, including those from space and SAR (search and rescue) radars [10,11,12,13,14,15,16]. However, these data are often unavailable for the specified time period, and thus numerical models and cutting-edge monitoring methods must be integrated. Moreover, in [1], it is stated that even the satellite-based systems can have a lot of false and missing detections. Out of the almost 800 synthetic aperture radar (SAR) images of spills observed in the Baltic Sea, approximately 94% were classified as having medium or low confidence levels. In contrast, only approximately 45% of the almost 70 confirmed incidents of oil pollution in the Latvian Baltic region were detected by SAR. Additionally, around 40% of the SAR detections corresponded to slicks or minor quantities of oil. By appropriately selecting a model that considers the influence of changing hydro-meteorological conditions on the movement and dispersion of pollution spots, relevant authorities can effectively respond with timely rescue actions in real-time. Such a model would enable accurate predictions and monitoring of the spread of pollutants, allowing for proactive measures to mitigate the impact on the environment and ensure swift response and containment efforts [2,9,11,15,17,18,19].

There are many methods allowing for the determination of the oil spill area, marking the domain, and predicting its trajectory to prevent ecological disasters. Their level of complexity, ease of use, and applicability vary depending on what features are most reasonable in a given sea basin. The big marine disasters concerned with oil spill [20] made the scientific community interested in this topic, and this resulted in the appearance of a large number of models [16,21,22,23]. According to the TAC (Technical Advisory Committee) on oil spill prevention and response it is possible to determine two categories of oil spill response [16,22,23,24]. Oil weathering models evaluate the changes over time without predicting the expected drift direction. The common approaches proposed in the literature, including the hydro-meteorological impacts, are often based on the historical statistical data, which is a reasonable approach. In the second group, there can be found models, such as: deterministic trajectory models, hind cast, stochastic models, and three-dimensional models (longitude, latitude, depth), which predict the oil behaviour and spread over time. Nevertheless, the chemical properties of oil were initially skipped [25] and later considered [26]. National Oceanic and Atmospheric Administration—NOAA provided a report, also including hazmat mapping as a third oil spill modelling tool [24]. A substantial amount of literature and numerous existing models have been developed to describe the movement of oil on the water surface, employing varying degrees of physics-based approaches. Examples of commonly utilized models include OpenDrift, SpillCalc, OSCAR, Mike, and SIMAP [1]. These models have been extensively studied and documented, providing valuable insights into understanding and predicting the behaviour of oil discharges in marine environments. The mathematical models and stochastic simulation software, such as GNOME, PISCES II, MEDSLINK-II, POSEIDON, STW, MOTHY, DieCAST-SSBOM, and COZOIL, can be used to mitigate the consequences of oil spills [6]. Several investigations have been performed on this topic [27,28,29,30,31,32,33,34,35]. Not all the studies are considered in the references because the research topic is extensive, but they are mostly based on ready-made models, instead of proposing the original ones from scratch.

An important issue in the implementation of the action of removing and stopping the spill is to determine precisely the shape and direction of oil spread [1,2,6,8,11,33,35,36,37]. To do this, author proposed their own written software based on modified and developed mathematical models [18,38,39,40] for predicting the horizontal spill movement and spreading in changing hydro-meteorological conditions through determining the oil spill domain composed of changing in time its surface elliptical sub-domains. Thus, the simple model of a surface oil slick as an ellipse is driven in the paper by a time-series of selected hydro-meteorological factors that change at random times. The algorithm developed for predicting oil spill trajectories and investigating the movement of domains is applicable universally, as it can be utilized in various areas of the Baltic Sea that share similar conditions related to hydrology and meteorology. This universality can be achieved through scientific experimentation and statistical identification of the model’s required parameters. By employing a simulation-based approach, real-time information can be provided to oil spill responders, enabling them to take prompt and informed actions. The research results obtained are the initial results. The model can be developed to provide a helpful reflection of reality with greater precision.

2. Process of Changing Hydro-Meteorological Conditions

2.1. Definition

Using the semi-Markov model [41,42], to mathematically model the process A(t), t ∈ 〈0, T〉, of varying hydro-meteorological (h-m) conditions during time T, T > 0, the necessary assumptions and denotations are presented in Figure 3 [43,44,45,46].

The paper considers the semi_-Markov approach as a more flexible and adaptable alternative to the conventional Markov model. Unlike the Markov model, which assumes exponential distributions for conditional sojourn times in specific h-m states, the semi-Markov approach allows for non-exponential distributions. This means that the model can accommodate any distribution of sojourn times in the respective states, providing a more realistic and accurate representation of the process. As a result, the semi-Markov approach is considered more practical, offering a better depiction of the real-world conditions.

In the following Section 2.2 and Section 2.3, the parameters for selected Baltic Sea waterway area are identified and the input data for Monte Carlo simulation are presented, assuming the initial time t = 0 (starting point). Once statistical data are collected, such as the number of initial states, counts of transitions between the states, and realizations of conditional sojourn times, the next step is to identify statistical parameters. This involves calculating the related probabilities and verifying the hypotheses regarding the distribution functions of the considered variables.

2.2. Parameters Identification for Selected Baltic Sea Waterway Area

Following a maritime accident and the subsequent release of oil, the resulting oil slick disperses and spreads across the water surface [7,20,22,26,35,47,48,49]. This happens due to gravity and interfacial tension [1,26,33,35,39]. The gravity influence decreases with time, but nevertheless the oil is still spreading [50]. Moreover, the hydro-meteorological conditions often significantly accelerate and lengthen the oil spread on the surface dependently on the wind direction and the wave height [1,23,26,33,49,51,52], especially in the open sea area. After consulting with experts from the polish Institute of Meteorology and Water Management, National Research Institute—IMWM-NRI, certain h-m factors were identified as having a significant impact on the trajectory of spill in the Baltic Sea waterway area. These selected factors were considered essential for understanding and accurately predicting the behaviour of the spills in that specific region:

• the velocity and direction of the wind;
• the height of the sea water level and waves;
• the direction of the currents;
• any obstacles to visibility such as fog or icing.

Therefore, we are engaged in studying the impact of these certain parameters on the shape and horizontal movement of the oil spill. We distinguish m = 6 states from {1, 2, …, 6} of the process of changing h-m conditions A(t), where t ∈ 〈0, T〉 [40,42] for the selected Baltic Sea open waters area defined by two parameters, the wave height (wh) and the wind speed (ws), presented in

Table 1. States of the process of varying h-m conditions for selected Baltic Sea open waters area.

States	wh [m]	ws [m/s]
1	0–2	0–17
2	2–5	0–17
3	5–14	0–17
4	0–2	17–33
5	2–5	17–33
6	5–14	17–33

.

Statistical methods are essential for determining the unknown parameters of the process to apply practically the simulation to the prediction of the oil discharge trajectory, considering the states of this process [40,42,53]. If the statistical data come from various experiments or the datasets are collected at different times but contain measurement points close to each other at the considered area, before the identification of parameters, the investigation of these data uniformity is necessary [42,54]. The uniformity testing procedure for the process of varying h-m conditions data, including realisation from different experiments or the datasets, is presented, e.g., in [42]. In the aforementioned work, the methods for estimating the probabilities of the initial process states (from Table 1), the probabilities of transitions between these states, and the distributions of the process’s sojourn times are discussed [17,18,40,41,42,44,45,46,47,48] and should be firstly applied.

The study focuses on analyzing the variations in h-m conditions specifically in the selected Baltic Sea waterway within the open waters area (excluding ports area and restricted water area) presented in Figure 4. The blue line is a waterway of a passenger ro-ro ship sailing from Gdynia to Karlskrona. The four points are the measurement points given from IMWM-NRI. The approximate coordinates are given on the map, whereas the four different measurement points: 1353, 1389, 1422, and 1458, are marked precisely in cooperation with the Institute.

Datasets of statistical nature were gathered over a six-year period (Θ = 6), specifically during the month of March, at open waters area [56]. First of all, the IMWM-NRI prepared data, including wh and ws (see Table 1) at the selected measurement points (1–4), as shown in Figure 4. The interval was 3 h. Thus, each year there were approximately 365 × 8 = 2920 records saved. Multiplying this number accordingly by 6 years was sufficient to perform the analysis. It made sense to choose one month due to varying conditions in the Baltic Sea, e.g., in June there are different winds than in February and March. This month was selected because from January to March there are more accidents than in the summer months, which is closely related to the storms, temperature, humidity, and visibility that occur. During that period, the weather in Sweden undergoes rapid changes, ranging from strong winds and storms to calm breezes. This factor is of significant importance in the investigation. The data were obtained from the IMWM-NRI under a copyright provision such that they can be presented in a processed form, not the exact original. Following the uniformity testing, the data collected from the four points (as depicted in Figure 4) were combined and subjected to analysis using methods outlined in [40,41,42,57,58,59,60,61]. Based on this consolidated dataset, the unknown fundamental parameters of the semi-Markov model of the process at the selected water area were assessed and evaluated.

Every day at 3 am in March during the 6-year period, the wh and the ws were checked. Then, according to Table 1, the initial state of the process, from {1, 2, …, 6}, was assigned. After that, all the initial realizations n_i, i = 1, 2, …, 6, were counted and presented in

Table 2. Numbers of initial realisations of the process.

State	1	2	3	4	5	6	Total
ni	405	237	0	0	27	11	680

.

Next, there were calculated the initial probabilities of the considered h-m conditions:

p₁(0) = 405/680 ≈ 0.595, p₂(0) = 0.349, p₃(0) = p₄(0) = 0, p₅(0) = 0.04, p₆(0) = 0.016;

(1)

Further, all the possible transitions between the 6 states were calculated. Every 3 h (3:00, 6:00, 9:00, 12:00, 15:00, 18:00, 21:00, 24:00), an evaluation was performed to assess whether the process continued or changed its state, i.e., if the two considered parameters (wh and ws) remain the same as 3 h ago or their value changed noticeably.

Table 3. The numbers of transitions between the states.

State	1	2	3	4	5	6
1→	–	1516	0	1	31	0
2→	1001	–	0	0	435	0
3→	0	11	–	0	1	0
4→	0	0	0	–	1	0
5→	0	298	2	0	–	88
6→	0	45	16	0	48	–

illustrates the obtained results.

According to the above table, there were calculated the probabilities of transitions between the process’ states:

p₁₂ = 1516/(1516 + 1 + 31) ≈ 0.98, p₁₅ = 0.02, p₂₁ = 0.7, p₂₅ = 0.3, p₃₂ = 0.92, p₃₅ = 0.08, p₄₅ = 1,
p₅₂ = 0.76, p₅₃ = 0.01, p₅₆ = 0.23, p₆₂ = 0.41, p₆₃ = 0.15, p₆₅ = 0.44.

(2)

The process is analyzed with respect to changes in time states. According to Section 2.1, the symbol θ_{kiki + 1} denotes the process’ random conditional sojourn time at the state k_i under the condition that its next transition will be to the state k_{i + 1}, i = 1, 2, …, n. For instance, the process is at the state k₁ = 1 (wh up to 2 m and small ws up to 17 m/s), but under the condition it will go to the state k₂ = 4 (the wave remains the same, but wind is changing). That change has happened once and the realisation time t_k1k2 = t₁₄ = 6 h. Next, we are at the state k₂ = 4 and the next transition was to the state k₃ = 5 (wave height increased up to 5 m plus there is a strong wind) and lasted for 3 h, which means t_k2k3 = t₄₅ = 3 h. An example of a process’s A(t) realization is presented in Figure 5.

Moreover, the hypotheses on the distribution functions W_kiki+1(t) were verified based on the sufficiently large realization number t_kiki+1, and the empirical distribution functions were determined for the numbers of realizations less than 30.

The collection of realizations of θ_kiki+1 representing the conditional sojourn times, where {k_i, k_i+1} = {{1,2}, {1,5}, {2,1}, {2,5}, {5,2}, {5,6}, {6,2}, {6,5}} were significantly large, and after verification it was confirmed that they exhibit chimney distributions (θ₁₂, θ₁₅, θ₂₁, θ₂₅, θ₅₂, θ₅₆, θ₆₂) and gamma distribution (θ₆₅), while the distribution functions of θ₁₄, θ₃₂, θ₃₅, θ₄₅, θ₅₃, θ₆₃ have the empirical distribution functions (Figure 6 and Figure 7). Formulae of the functions W_{kiki + 1}(t) are approximated to three decimal places, whereas in the graphs there are illustrated functions drawn without that approximation. The chimney distribution is introduced in a book [42], Chapter 2, page 58. This distribution got its name from chimneys because of the chimney look of the density function. The chimney distributions were used because they fit very well to the given data sets. However, a few distribution functions have the empirical distribution functions and none of the common distributions can fit well to the data (at least 30 distributions were checked with 1% to 10% critical interval). Moreover, some transitions to next states do not appear frequently and there were a few realisations. For instance, a transition between states 1 → 4: wh up to 2 m and small ws up to 17 m/s that increases above that value without much change in the wave height happened once during the considered period of time. Thus, empirical distributions were used.

The illustrated graphs and the presented formulae for distribution functions in Figure 6 and Figure 7 are later used in Section 2.3 to predict the h-m conditions using simulation method. There are multiple methods available for generating random samples from a specified probability distribution. Examples of such transform methods include the inverse method, the Box-Muller (B-M) method, and Marsaglia and Tsang’s rejection sampling (MTRS) method. This paper employs the inverse transform method for generating the realizations from $t_{k_i{k}_{i+1}}^{(i)}(h):=W_{k_i{k}_{i+1}}^{-1}(h)$. The gamma distribution realisations are generated using the MTRS method. The values are generated from transformation of the Gaussian function.

2.3. Input Data for Simulation

The data input to the simulation method based on the parameters identified in Section 2.2 are as follows:

• the experiment time T = 48 h;
• formula for generating the initial state coming from (1)

$$k_i\left(q\right)=\{\begin{array}{ll}1,\hfill & 0\le q<0.595\hfill \\ 2,\hfill & 0.595\le q<0.944\hfill \\ 5,\hfill & 0.944\le q<0.984\hfill \\ 6,\hfill & 0.984\le q<1,\hfill \end{array}$$

(3)

where the state is selected by generating from the distribution of the initial probabilities determined by (1) and q is a number that is randomly generated from the uniform distribution over the interval 〈0, 1);

• formula for generating next states

$$k_{i+1}\left(g\right)=\{\begin{array}{l}\{\begin{array}{l}2,\hspace{1em}0\le g<0.98\hfill \\ 5,\hspace{1em}0.98\le g<1\hfill \end{array}\hspace{1em}{k}_i=1\hfill \\ \{\begin{array}{l}1,\hspace{1em}0\le g<0.7\hfill \\ 5,\hspace{1em}0.7\le g<1\hfill \end{array}\hspace{1em}{k}_i=2\hfill \\ \{\begin{array}{l}2,\hspace{1em}0\le g<0.92\hfill \\ 5,\hspace{1em}0.92\le g<1\hfill \end{array}\hspace{1em}{k}_i=3\hfill \\ 5,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{k}_i=4\hfill \\ \{\begin{array}{l}2,\hspace{1em}0\le g\le 0.76\hfill \\ 3,\hspace{1em}0.76\le g\le 0.77\hfill \\ 6,\hspace{1em}0.77\le g<1\hfill \end{array}\hspace{1em}{k}_i=5\hfill \\ \{\begin{array}{l}2,\hspace{1em}0\le g\le 0.41\hfill \\ 3,\hspace{1em}0.41\le g\le 0.56\hfill \\ 6,\hspace{1em}0.56\le g<1\hfill \end{array}\hspace{1em}{k}_i=6,\hfill \end{array}$$

(4)

constructed using (2), where g is a number that is randomly generated from the uniform distribution over the interval 〈0,1);

• formulae to enable determination of realizations of the empirical conditional sojourn times of the process that exhibits varying h-m conditions are presented in Figure 6 and Figure 7.

The experiment time is chosen as a short-term forecast, because the first two days of an oil discharge are of utmost importance to accurately determine the shape and direction of spread to implement quickly the action of removing the pollution. The simulation time can of course be increased, even by much. However, for long time prediction, another method should be used.

3. Algorithm

Figure 8 presents an expanded algorithm that proposes a method for determining the trajectory of oil discharge by taking into account its horizontal movement, which is influenced by varying h-m conditions. This is based on a probabilistic approach [17,18,40] and an extended approach from [39].

The process of varying h-m conditions at the selected water location where the oil spill occurred, depending on the specific incident being referred to, is constructed in the previous section. The h-m conditions presented in Table 1 are considered in the input model in Equations (3) and (4) which are generated from Equations (1) and (2) and equations derived from Figure 6 and Figure 7. Moreover, Equation (5) determines the parametric equations for the drift trend curve of the central point of the oil spill domain under different h-m conditions. These equations serve as the input for the model, allowing for the prediction of how oil drifts over time based on the specific conditions experienced.

The algorithm can consider more than two hydro-meteorological factors, e.g., wave height, wind speed, wind direction, ice conditions, sea water level, etc., and generates surface elliptical sub-domains, where actions of mitigation should be taken to reduce the pollution at the oil spill water area. In this paper, only two parameters were intentionally taken into consideration. For example, to add a third, fourth, or more parameters to Table 1, the data should be firstly grouped into sets and the table needed to be extended accordingly. For the purpose of this article, in order not to present (already) too many mathematical calculations, a simple example is illustrated in Section 4, to make the paper familiar to the reader.

This algorithm operates under the assumption that the oil spill central point $(X^{k_i}(t),Y^{k_i}(t)),$ at t, t ∈ 〈0, T〉, while the process of varying h-m is at the state k_i, i = 1, 2, …, n, and has the normal distribution with the parameters that will be slightly modified [39,40]:

• expected values $m_X^{k_i}(t),$ $m_Y^{k_i}(t);$
• standard deviations ${\sigma }_X^{k_i}(t)$, ${\sigma }_Y^{k_i}(t)$;
• correlation coefficient ${\rho }_{XY}^{k_i}(t)$.

Moreover, it is assumed that the oil spill domain radius $r^{k_i}(t)$ at k_i, i = 1, 2, …, n, is dependent on time t ∈ 〈0,T〉.

The points ($m_X^{k_i}(t),$ $m_Y^{k_i}(t)$), t ∈ 〈0,T〉, k_i = 1, 2, …,~m, create a curve K^k called an oil spill central point drift trend (CPDT)

$$K^{ki}:\{\begin{array}{l}{x}_{}^{k_i}=x_{}^{k_i}(t)\\ y_{}^{k_i}=y_{}^{k_i}(t),t\in <0,T>,\end{array}$$

(5)

where x, y are given in meters. An exemplary CPDT is presented in Figure 9.

Therefore, assuming CPDT at each state the sequence of oil spill domains is given as follows

$${\overline{D}}^{k_i}((s_{i-1}+1)\Delta t),{\overline{D}}^{k_i}((s_{i-1}+2)\Delta t),\dots ,{\overline{D}}^{k_i}(s_i\Delta t),$$

(6)

with the step of time ∆t, where ${\overline{D}}^{k_i}(t)$, for t equals to (s_i–1 + 1)∆t, (s_i–1+ 2)∆t, …, s_i∆t, are given by

$$\begin{array}{l}{\overline{D}}^{k_i}(t)=\{(x,y):\frac{1}{1-{({\rho }_{XY}^{k_i}(t))}^2}[\frac{{(x-m_X^{k_i}(t))}^2}{{({\sigma }_X^{k_1}(t))}^2}–2{\rho }_{XY}^{k_i}(t)\frac{(x-m_X^{k_i}(t))(y-m_Y^{k_i}(t))}{{\sigma }_X^{k_i}(t){\sigma }_Y^{k_i}(t)}\\ +\frac{{(y-m_Y^{k_i}(t))}^2}{{({\sigma }_Y^{k_i}(t))}^2}]\le c^2\},t\in \langle 0,T\rangle ,k_i=1,2, \dots ,m,i=1,2, \dots ,n.\end{array}$$

(7)

The sequence (6) is determined for the experiment time T in the successive intervals ((s_i–1 + 1)∆t, s_i∆t〉 at state k_i, with s_i that satisfy the inequality

$$\left(s_i{}_{-1}\right)\Delta t<{\displaystyle \sum _{j=1}^i}{t}_{kj kj+1}\le s_i\Delta t,i=1,2, \dots ,n,$$

(8)

where $s_n\Delta t\le T,$ and t_{kj kj + 1} = $W_{k_j{k}_{j+1}}^{-1}$(h), j, = 1, 2, …,i, i = 1, 2, …, n, are the inverse functions of the conditional distribution functions W_{kj kj + 1}(t) and h is a number that is randomly generated over the interval 〈0, 1).

In the experiment, the considered domain is described by combining the determined domains of the sequences (6). The resulting domain represents the total extent under investigation:

$${\overline{\overline{\mathit{D}}}}^{k_1,k_2,\dots ,k_n}={\displaystyle \underset{i=1}{\overset{n}{\cup }}{\displaystyle \underset{j=1}{\overset{s_i}{\cup }}{\overline{D}}^{k_i}\left((s_{i-1}+j)\Delta t\right)}},$$

(9)

with:

• c² = −2ln(1 − p);
• expected values

$$m_X^{k_i}(t):=m_X^{k_{i-1}}(s_{i-1}\Delta \mathrm{t})+m_X^{k_i}(a_i\Delta \mathrm{t}),$$

(10)

$$m_Y^{k_i}(t):=m_Y^{k_{i-1}}(s_{i-1}\Delta \mathrm{t})+m_Y^{k_i}(a_i\Delta \mathrm{t});$$

(11)

• standard deviations

$${\overline{\sigma }}_X^{k_i}(t):={\sigma }_X^{k_i}((s_{i-1}+a_i)\Delta \mathrm{t})+{\displaystyle \sum _{j=1}^i{r}^{k_j}}(b_j\Delta \mathrm{t}),$$

(12)

$${\overline{\sigma }}_Y^{k_i}(t):={\sigma }_Y^{k_i}((s_{i-1}+a_i)\Delta \mathrm{t})+{\displaystyle \sum _{j=1}^i{r}^{k_j}}(b_j\Delta \mathrm{t}),$$

(13)

with radiuses

$$r^{k_j}(\mathrm{t}):=r^{k_j}(b_j\Delta \mathrm{t}),j=1,2, \dots ,i;$$

(14)

• correlation coefficients ${\rho }_{XY}^{k_i}(t)$,

for a_i = 1, 2, …, b_i, b_i = 1, 2, …, s_i − s_i–1, i = 1, 2, …, n.

4. Application

Based on the statistical data presented in Section 2.2 and utilizing the input data from Section 2.3, the application of the algorithm outlined in Section 3 becomes possible. By employing the Monte Carlo approach, as described in references [19,26,35,62,63,64], it becomes feasible to predict the changes in the shape of the oil spill domain and forecast the trajectory of the oil under varying h-m conditions for the open water area of the Baltic Sea.

4.1. Experiment

We arbitrarily assume that the points $(m_X^{k_i}(t),m_Y^{k_i}(t))$, t ∈ 〈0,48〉, k_i ∈ {1, 2, …, 6}, i = 1, 2, …, n, create an oil spill CPDT composed of curves K^ki, i = 1, 2, …, n [40,65,66]. According to (5), we have

$$K^{ki}:\{\begin{array}{l}{x}_{}^{k_i}=t^{k_i}\\ y_{}^{k_i}=t,t\in ((s_{i-1}+1)\Delta t,s_i\Delta t\rangle ,k_i\in \{1,2,\dots ,6\},i=1,2,\dots ,n.\end{array}$$

(15)

By repeatedly applying the procedure outlined in Section 3 with a time step of ∆t = 1 h, the domain that we are looking for can be obtained. This can be determined based on arbitrarily assumed standard deviations and the correlation coefficient of the coordinates of the oil spill CPDT and the radiuses of the domain, as specified in references [50,67,68]:

• ${\sigma }_X^{k_i}(t)$ = ${\sigma }_Y^{k_i}(t)$ = ${\sigma }_{}^{k_i}(t)$ = 0.1 + 0.2∙t;
• ${\rho }_{XY}^{k_i}(t)$ = 0.8;
• $r^{k_i}(t)=$0.5 + 0.5∙t,

for t ∈ 〈0,48〉, k_i ∈ {1, 2, …, 6}, i = 1, 2, …, n, that in the real practice should be statistically identified using the methods given in [40].

First, a random number q is generated from the uniform distribution on the interval 〈0,1), and in this case q ≅ 0.62. Next, the initial state k₁, where k₁ ∈ {1, 2, …, 6}, is selected according to (3), resulting in k₁(0.62) = 2.

Subsequently, another random number g is drawn, and in this case, g ≅ 0.45. With the fixed k₁ = 2, we determine the next state k₂, where k₂ ∈ {1, 2, …, 6}, using Equation (4), i.e., k₂(0.45) = 1.

Afterwards, a random number h is drawn, and in this case, h ≅ 0.20. Considering the fixed states k₁ = 2 and k₂ = 1, we generate the first realization t_{kiki + 1} = $t_{k_1{k}_2}^{(1)}$ of θ_k1k2 from an assigned probability distribution, as indicated in Figure 6 and Figure 7. In this example, $t_{k_1{k}_2}^{(1)}$ = t₂₁ = 6.78. Consequently, we have (s₁ − 1) < 6.78 ≤ s₁.

Hence, s₁ = 7 and s₀ = 0, s₁ − s₀ = s₁ = 7. Subsequently, we compare the value of s₁ with time T, which is equal to 48. Observing that s₁ = 7, which is less than 48, we proceed to draw the sequence of domains accordingly:

$$\begin{array}{l}{\overline{\overline{\mathit{D}}}}^{k_1}(b_1)={\overline{\overline{D}}}^{k_1}(s_0+a_1):={\overline{D}}^k(0+a_1)\\ =\{(x,y):\frac{1}{1-{(0.8)}^2}[\frac{{(x-m_X^k(a_1))}^2}{{({\overline{\sigma }}_{}^k(a_1))}^2}–2\cdot 0.8\frac{(x-m_X^k(a_1))(y-m_Y^k(a_1))}{{\overline{\sigma }}_{}^k(a_1)\cdot {\overline{\sigma }}_{}^k(a_1)}\\ +\frac{{(y-m_Y^k(a_1))}^2}{{({\overline{\sigma }}_{}^k(a_1))}^2}] \le -2ln(1-0.95)\cong 5.99\},\end{array}$$

(16)

where

$$m_X^k(a_1):=m_X^{k_0}(s_0)+m_X^{k_1}(a_1)={(a_1)}^2,$$

(17)

$$m_Y^k(a_1):=m_Y^{k_0}(s_0)+m_Y^{k_1}(a_1)=a_1,$$

(18)

$$\begin{array}{l}{\overline{\sigma }}_{}^k(a_1):={\overline{\overline{\sigma }}}_{}^{k_1}(s_0+a_1)={\sigma }_{}^{k_1}(s_0+a_1)+{\displaystyle \sum _{j=1}^1{r}^{k_j}}(b_j)={\sigma }_{}^{k_1}(a_1)+r^{k_1}(b_1)\\ =0.1+0.2a_1+0.5+0.5b_1=0.2a_1+0.5b_1+0.6,\end{array}$$

(19)

for a₁ = 1, 2, …, b₁, b₁ = 1, 2, …, 7, composed of the elliptical sub-domains shown in Figure 10.

Further, the substitution of k₂ = 1 is proceeded by randomly generating another set of numbers, g ≅ 0.76 and h ≅ 0.30. We select state k₃(0.76) = 2, and generate another realization, $t_{k_2{k}_3}^{(2)}$ = t₁₂ = 18.75, for the conditional sojourn time. Once we have these realizations, we calculate their sum:

$$t_{k_1{k}_2}^{(1)}+ t_{k_2{k}_3}^{(2)}=t_{21}+t_{12} \cong 6.78+18.75=25.53.$$

(20)

Considering the inequality (s₂ − 1) < 25.53 ≤ s₂, we deduce that s₂ = 26 and the difference s₂ − s₁ = 26 − 7 = 19. Next, we compare s₂ with time T and observe that s₂ = 26, which is less than 48. Consequently, we obtain the sequence:

$$\begin{array}{l}{\overline{\overline{\mathit{D}}}}^{k_1,k_2}(b_2)={\overline{\overline{D}}}^{k_2}(s_1+a_2):={\overline{D}}^k(7+a_2)\\ =\{(x,y):\frac{1}{1-{(0.8)}^2}[\frac{{(x-m_X^k(7+a_2))}^2}{{({\overline{\sigma }}_{}^k(7+a_2))}^2}\\ –2\cdot 0.8\frac{(x-m_X^k(7+a_2))(y-m_Y^k(7+a_2))}{{\overline{\sigma }}_{}^k(7+a_2)\cdot {\overline{\sigma }}_{}^k(7+a_2)}\\ +\frac{{(y-m_Y^k(7+a_2))}^2}{{({\overline{\sigma }}_{}^k(7+a_2))}^2}]\le 5.99\},\end{array}$$

(21)

of domains, where

$$m_X^k(7+a_2):=m_X^{k_1}(s_1)+m_X^{k_2}(a_2)=7^2+a_2,$$

(22)

$$m_Y^k(7+a_2):=m_Y^{k_1}(s_1)+m_Y^{k_2}(a_2)=7+a_2,$$

(23)

$$\begin{array}{l}{\overline{\sigma }}_{}^k(7+a_2):={\overline{\overline{\sigma }}}_{}^{k_2}(s_1+a_2)={\sigma }_{}^{k_2}(7+a_2)+r^{k_1}(b_1)+r^{k_2}(b_2)\\ =0.1+0.2\left(7+a_2\right)+0.5+0.5b_1+0.5+0.5b_2\\ =0.2a_2+0.5b_1+0.5b_2+2.5.\end{array}$$

(24)

for a₂ = 1, 2, …, b₂, b₂ = 1, 2, …, 19.

The sequence of oil spill domains in the time interval 〈s₁, s₂) = (7, 26〉 is partly illustrated for selected t = 8 h, 9 h, 15 h, 26 h, as shown in Figure 11. We can notice that it consists of elliptical sub-domains with bigger radiuses than those shown in Figure 10.

In the next step, we generate random numbers, g ≅ 0.03 and h ≅ 0.80. By selecting the state k₄(0.03) = 1, we then obtain a realization as $t_{k_3{k}_4}^{(3)}$ = 52.3. The entire sojourn time is:

$$t_{k_1{k}_2}^{(1)}+ t_{k_2{k}_3}^{(2)}+ t_{k_3{k}_4}^{(3)}=t_{21}+t_{12}+t_{21} \cong 6.78+18.75+52.3=77.83,$$

(25)

and (s₃ − 1) < 77.83 ≤ s₃.

As s₃ is greater than the experimental time 48 h, we substitute and calculate the difference s₃ − s₂ = 48 − 26 = 22. Based on this, we accordingly proceed to draw the sequence

$$\begin{array}{l}{\overline{\overline{\mathit{D}}}}^{k_1,k_2,k_3}(b_3)={\overline{\overline{D}}}^{k_3}(s_2+a_3):={\overline{D}}^k(26+a_3)\\ =\{(x,y):\frac{1}{1-{(0.8)}^2}[\frac{{(x-m_X^k(26+a_3))}^2}{{({\overline{\sigma }}_{}^k(26+a_3))}^2}\\ –2\cdot 0.8\frac{(x-m_X^k(26+a_3))(y-m_Y^k(26+a_3))}{{\overline{\sigma }}_{}^k(26+a_3)\cdot {\overline{\sigma }}_{}^k(26+a_3)}\\ +\frac{{(y-m_Y^k(26+a_3))}^2}{{({\overline{\sigma }}_{}^k(26+a_3))}^2}] \le 5.99\},\end{array}$$

(26)

of domains, where

$$m_X^k(26+a_3):=m_X^{k_2}(s_2)+m_X^{k_3}(a_3)=68+{(a_3)}^2,$$

(27)

$$m_Y^k(26+a_3):=m_Y^{k_2}(s_2)+m_Y^{k3}(a_3)=26+a_3,$$

(28)

$$\begin{array}{l}{\overline{\sigma }}_{}^k(26+a_3):={\overline{\overline{\sigma }}}_{}^{k_3}(s_2+a_3)\\ ={\sigma }_{}^{k_3}(26+a_3)+r^{k_1}(b_1)+r^{k_2}(b_2)+r^{k3}(b_3)\\ =0.1+0.2\left(26+a_3\right)+0.5+0.5b_1+0.5+0.5b_2+0.5+0.5b_3\\ =0.2a_3+0.5b_1+0.5b_2+0.5b_3+6.8,\end{array}$$

(29)

for a₃ = 1, 2, …, b₃, b₃ = 1, 2, …, 22.

The sequence of oil spill domains (composed of the elliptical sub-domains) in the time interval 〈s₃, s₄) = (26, 48〉 is partly illustrated for selected t = 27 h, 28 h, 38 h, 47 h in Figure 12. Figure 13 displays the final domain at t = 48 h. The sequence for varying h-m conditions after two days of the oil discharge is presented. The simulation time can be increased and for this purpose the algorithm should be continued, i.e., the process involves substituting i with j and subsequently drawing new random numbers g, h. This is conducted to select the states k_i+1 and generate different realizations $t_{k_i{k}_{i+1}}^{(i)}(h)$. This is iteratively repeated until the cumulative sum ${\displaystyle \sum _{j=1}^i}$t_kjkj+1 obtained from all the generated realizations $t_{k_i{k}_{i+1}}^{(i)}(h)$ reach a new increased experiment time T. Then, the necessary parameters can be calculated and by following the aforementioned iteration process, sequences of domains under varying h-m conditions can be obtained. In the experiment, the oil spill domain is represented by the cumulative sum of sub-domains of an ellipse shape.

4.2. Discussion of Results and Further Developments

The preliminary results are obtained for arbitrarily assumed expressions of expected values $m_X^{k_i}(t)$ and $m_Y^{k_i}(t),$ standard deviations ${\sigma }_X^{k_i}(t)$ and ${\sigma }_Y^{k_i}(t)$, correlation coefficients ${\rho }_{XY}^{k_i}(t)$ of the oil spill CPDT curves K^ki, its equation, and the radiuses $r_{}^{k_i}(t)$ of the moving oil spill elliptical sub-domains. In the real practice, the unknown parameters should be statistically identified.

Figure 13 presents the final domain for a single simulation trial. To enhance the algorithm and simulation procedure, it is feasible to repeat the trials multiple times and calculate the average domain. This approach allows for a more comprehensive analysis and a better understanding of the overall behaviour of the oil spill.

The author aims to conduct a practical scientific experiment in specific water areas of the Baltic Sea with the intention of statistically identifying the unknown parameters associated with the assumed elements of the model proposed. This ambitious endeavour seeks to determine the forms of these model components under fixed hydro-meteorological conditions in the sea area. The insights gained from this study can be used to predict spill domains in other sea water areas with similar conditions. However, the experiments at the considered sea water areas will consume much time and cost, especially if performed for different kinds of spills and also including the physical oil characteristics and other processes, such us evaporation, emulsification, and thickness of the oil spill layer. The simulation results of the oil slick domain in terms of moving trajectories can be developed to ensure the greatest precision and correspondence to reality as possible.

5. Conclusions

This paper reports work that is of interest for practical safety in wide water areas. It proposes the utilization of the simulation technique, which employs a stochastic approach, for determining the trajectory of oil spills and investigating their horizontal domain movement. A probabilistic model is introduced and identified for the process of changing hydro-meteorological conditions at the water surface. A two-dimensional stochastic process is then employed to determine the central point position of the spill domain, and parametric equations for its trend curves under different h-m conditions are defined. Using the aforementioned approach, a comprehensive simulation algorithm is formulated to simulate the stochastic movement of oil spills within a dynamic hydro-meteorological environment. This algorithm takes into account the changing conditions over time. The practical application is demonstrated by applying it to the open water area of the Baltic Sea as a specific case study. The study specifically emphasizes two significant parameters: wave height and wind speed. By considering these factors, the algorithm provides insights into the movement patterns of oil in response to variations in h-m conditions. After successful uniformity testing and joining realisations of conditional sojourn times from different experiments or datasets into new sets of data, the unknown parameters investigation methods were applied and the simulation technique used. The prediction of oil spill movement and spreading in varying hydro-meteorological conditions is performed through determining the trajectory and movement of the oil spill domain composed of changing its surface elliptical sub-domains in time. The final result is obtained after reaching a fixed experiment time T.

The model presented in this paper serves as a valuable tool in the case of spills in areas of various oceans, seas, and inland water. The proposed simulation algorithm for oil spill domain movement investigation is universal in the sense that we can apply it in similar areas of the Baltic Sea with comparable hydro-meteorological conditions after the necessary identification of statistical model parameters through the scientific experiment. Transferring the findings to other water regions that share similar h-m conditions would undoubtedly enhance and optimize environmental protection efforts. This method provides a rapid and precise means of determining oil discharge trajectories and predicting their domain movement. The research outcomes obtained from this study serve as initial results in the field of oil spill simulation.

As oil spills can have a profound impact on the reliability of operations, profitability of investments, and operating costs associated with various industries and sectors [69,70], it seems reasonable to devote significant attention and resources to prevent and effectively respond to oil discharges, ensuring the sustainability and resilience of these industries in the face of potential environmental risks. Moreover, the aviation industry, with its reliance on fuel and oil can have the potential impact on oil spills [71,72]. Investing in early warning systems and improving spill response capabilities can help minimize the impact of e.g., flood-induced oil spills [73]. The method proposed in the paper can be applied to various types of spills, thereby mitigating the detrimental consequences for the environment.