Operational Ocean Modelling in Support of Forensic Investigations: A Backward Lagrangian Drift Modelling for Migrant Shipwreck Reconstruction

1. Introduction

Marine and oceanic currents may transport passive drifters (including marine plankton, buoys, dinghies, and floating debris) for long distances. For this reason, oceanographic and marine circulation models may be employed by governmental and humanitarian agencies, such as the INHI (Italian Navy Hydrographic Institute), MOAS (Migrant Offshore Aid Station), UNHCR (United Nations High Commissioner for Refugees), and IOM (International Organisation for Migration), to simulate the probable trajectories of drifting objects and, in some cases, to reconstruct their possible routes retrospectively.

Most studies concerning drifter trajectory prediction have focused on biological (particularly zooplankton and marine turtles) or environmental targets [1,2,3]. The application of such predictive models to human remains is comparatively recent. The first study in which a drifting object was represented by a human body was published only in 2002 [4]. Since then, additional investigations have highlighted the potential forensic relevance of oceanographic drift modelling, particularly for estimating the probable origin, trajectory, and recovery area of drowned individuals at sea [5]. These approaches may therefore provide valuable support in forensic investigations involving maritime disasters, missing persons, and migrant shipwrecks. In this last case, over the last few decades, irregular migration flows from North Africa and transcontinental Turkey toward the European countries bordering the Mediterranean Sea (particularly Italy, Spain, France, Malta, and Greece) have increased substantially [6,7]. According to the main governmental and humanitarian agencies, and mass media news, the principal known maritime routes across the Central Mediterranean are the most dangerous in the world [8] and include the corridors connecting Morocco and northwestern Algeria to Spain and the Balearic Islands; northeastern Algeria to Sardinia; Tunisia to Pantelleria and Lampedusa Island; Libya to southern Sicily, including Lampedusa and Pantelleria Islands and the southern Sicilian coastline; Egypt to southeastern Sicily, Calabria, and Crete; and Turkey to Calabria. Secondary migration routes also connect Tunisia with southern Sardinia. These crossings are typically undertaken using overcrowded (from 15–20 to ∼50 individuals, including women and children) and structurally inadequate vessels. Consequently, due to both the very dangerous maritime conditions and the increasing intensity of migratory movements, shipwrecks occur frequently, causing the deaths by drowning of thousands of migrants each year. As a result of this ongoing humanitarian crisis, the Mediterranean Sea has increasingly been described as a ”cemetery”, where human remains are often recovered accidentally by fishermen, sailors, or recreational boaters, sometimes months after a maritime disaster. These migration journeys, notwithstanding the existence of several international and transcontinental protocols, conventions, and humanitarian actions are still very frequent, often without a clearly defined final destination.

Within this complex geopolitical and forensic context, in March–April 2024, several human bodies were recovered in the Italian waters of the southern Tyrrhenian Sea, around the Aeolian Archipelago and the northeastern coast of Sicily near the Strait of Messina. The recoveries received considerable attention in the national and regional media, which reported the hypothesis that the victims might have originated from a migrant shipwreck that reportedly occurred on 29 March 2024 in the Sicilian Channel, offshore Lampedusa Island (Agrigento, southern Sicily [9]). Overall, between March and June 2024, five human bodies were recovered in the maritime sector of Messina.

However, considering the considerable distance between the presumed shipwreck location and the recovery sites, as well as the complexity of the Mediterranean circulation patterns, an important forensic and oceanographic question emerged: was it physically plausible for bodies originating from a shipwreck near Lampedusa to drift northward a considerable distance from the scene of death, in the Ionian Sea, along the eastern Sicilian coastline, and eventually reach the Tyrrhenian Sea, overtaking the Strait of Messina? Preliminary analyses based on Copernicus Marine Service data [10], including information on prevailing sea currents and wind fields in the central Mediterranean during the relevant period, appeared inconsistent with such a drift scenario. This discrepancy raised the possibility that the recovered individuals may have originated from a different migratory route or from another undocumented maritime accident.

Reconstructing the most probable drift trajectories of migrant shipwrecks may provide crucial investigative support to judicial authorities, forensic pathologists, and international security agencies. In particular, spatio-temporal predictive modelling can assist in narrowing down possible shipwreck locations, identifying migratory routes, supporting victim identification processes, and guiding intelligence and rescue investigations.

Within the framework of the present study, an assessment was carried out to investigate whether open-access data and publicly available tools can represent a reliable and effective support for the identification of the probable departure location of migrant vessels involved in shipwreck events. Particular attention was devoted to evaluating the capability of these resources to reconstruct the trajectories followed by drifting bodies and to relate the final recovery locations to potential source areas. The analysis was motivated by the fact that in many documented cases, the exact location of the maritime incident and the vessels’ original departure points remain unknown or only partially constrained. In this context, integrating publicly available oceanographic and meteorological datasets with numerical tools and trajectory analyses can help narrow the range of possible origin locations. The proposed approach aims not only to support forensic and humanitarian investigations on missing migrants at sea but also to contribute to a better understanding of maritime migration routes, drift dynamics in the central Mediterranean Sea, and the spatial distribution of shipwreck-related events. Furthermore, the availability of reproducible, openly accessible methodologies may facilitate future applications by researchers, operational agencies, and humanitarian organisations engaged in search, rescue, and identification activities.

The rest of the paper is organised as follows. Section 2 presents the adopted methodology, including the available data, the numerical model used to reconstruct the possible shipwreck location, the simulated scenarios and the identification of the eventual routes. Section 3 reports the obtained results, while the paper ends with discussion and conclusions.

2. Materials and Methods

The proposed framework combines oceanographic observations, atmospheric forcing, and Lagrangian numerical simulations to reconstruct plausible transport pathways and to identify potential departure regions associated with undocumented maritime accidents (Figure 1).

By coupling reverse drift simulations with shipwreck information, the methodology aims to reduce uncertainties in the reconstruction of the events leading to the discovery of bodies along coastlines or in offshore areas. Furthermore, the analysis of meteomarine data within the inferred source regions provides additional insight into the environmental factors that may have contributed to vessel instability, capsizing, or sinking events. A correct identification of the accident site could also help in locating additional victims or the shipwrecked vessel. Specifically, the methodology developed and implemented in the present study consists of the following main steps:

1.

identification and georeferencing of the body recovery locations, including the collection of available forensic data used to constrain the temporal window of the simulations;

2.

selection of the OpenDrift Leeway model as the numerical framework for the backwards-in-time simulation of drifting trajectories, on the basis of its suitability for simulating the drift of persons and objects at the sea surface under the combined effect of surface currents, wind, and Stokes drift;

3.

collection and pre-processing of high-resolution oceanographic and meteorological datasets from the Copernicus Marine Service (CMEMS), comprising sea surface currents, wind fields, and Stokes drift, used as environmental forcing for the Lagrangian simulations;

4.

backwards-in-time simulation of multiple possible release scenarios of virtual drifters starting from the recovery locations;

5.

spatial and temporal analysis of the simulated trajectories to identify the most probable origin areas and constrain the spatiotemporal window of the shipwreck event;

6.

cross-referencing of the high-probability regions with documented shipwreck incidents to identify the most plausible event associated with the recovered bodies;

7.

analysis of the wave climate along the estimated transit corridor of the vessel to assess whether severe sea state conditions consistent with the capsizing of the boat were present during the estimated time window of the incident.

2.1. Identification of the Body Recovery Locations

A total of five human remains were recovered in the area between the Aeolian Islands and the northeastern coast of Sicily between March and June 2024.

The first body (ID 1) was found on 17 March 2024 at approximate coordinates 38°19’52"N, 14°59’03"E, offshore Vulcano Island; the individual was of Caucasian appearance. The second body (ID 2) was recovered on 13 April 2024 at coordinates 38°16.607’N, 14°59.974’E, in the open sea between the Gulf of Patti and Vulcano Island; the individual was Caucasian, presented a tattoo, and was wearing dark cargo trousers, with barnacle (Lepas spp.) colonisation noted on the remains. The third body (ID 3) was recovered on 14 April 2024 at coordinates 38°33.2’N, 14°36.9’E, offshore Filicudi Island; the individual was Caucasian, presented a tattoo, and was wearing an orange life jacket, with barnacle colonisation also noted. The fourth body (ID 4) was found on 18 April 2024 at coordinates 38°16’14.52"N, 15°28’43.43"E, on the beach of Rodia; the individual was Caucasian and presented tattoos on the left thigh and left flank. The fifth body (ID 5) was recovered on 30 June 2024 at coordinates 38°33.873’N, 14°41.046’E, in the open sea between Salina and Filicudi islands; the individual was Caucasian and was found fully clothed and wearing a life jacket.

Based on the Lepas spp. colonisation observed on bodies ID 2 and ID 3, the forensic pathologist estimated a minimum time of subimmersion of approximately 2–3 months. The minimum subimmersion interval was also confirmed by the transformative cadaveric phenomena taking into account extrinsic factors (environmental and climatic) This forensic estimate was used to define the temporal window of the backwards simulations. A summary of recovery locations, dates, and available forensic data is provided in Table 1.

Figure 2 shows the study area and the five body recovery locations.

2.2. Numerical Model

To identify the region where the shipwreck occurred, numerical simulations were conducted using the OpenDrift model, developed by the Norwegian Meteorological Institute (MET Norway). OpenDrift is an open-source Python framework for modelling Lagrangian trajectories of particles, such as ocean drifters, floating objects, or substances suspended in the atmosphere or ocean. It is highly modular and flexible, enabling integration of data from ocean, atmospheric, and wave models, as well as direct measurements [11].

The Leeway module in OpenDrift is specifically designed to simulate the motion of wind-affected floating objects. It is based on operational models developed by the U.S. Coast Guard [12,13] and further adapted by MET Norway. The model has been used in several studies (e.g., [14,15]).

The module simulates the drift of floating objects by combining wind-induced leeway, ambient ocean currents, and surface wave-induced (Stokes) drift. The term leeway refers to the drift of a floating object caused by surface winds acting on the portion of the object above the waterline. The leeway velocity (${\overrightarrow{V}}_{\mathrm{leeway}}$) is decomposed into two components relative to the wind direction: a downwind component ($V_{dw}$), describing motion in the wind direction, and a crosswind component ($V_{cw}$), accounting for lateral deviation due to the object’s shape, balance, and hydrodynamic characteristics.

These components are determined from empirically derived coefficients specific to each object type, with added stochastic perturbations to account for natural variability:

$$V_{dw}=(S_{dw}+{ϵ}_{dw})U_{10}+O_{dw}$$

(1)

$$V_{cw}=(S_{cw}+{ϵ}_{cw})U_{10}+O_{cw}$$

(2)

where S represents the slope coefficients, O the offsets, $ϵ$ the stochastic perturbations, and $U_{10}$ the wind speed at 10 m above sea level.

The total velocity of a particle also includes the contribution of ambient surface currents ($V_{\mathrm{curr}}$) and Stokes drift ($V_{\mathrm{Stokes}}$), i.e. the net mass transport velocity of surface fluid parcels induced by the propagation of gravity waves, whose components are:

$${\displaystyle \frac{dX}{dt}}=V_x^{\mathrm{leeway}}+V_x^{\mathrm{curr}}+V_x^{\mathrm{Stokes}}$$

(3)

$${\displaystyle \frac{dY}{dt}}=V_y^{\mathrm{leeway}}+V_y^{\mathrm{curr}}+V_y^{\mathrm{Stokes}}$$

(4)

2.3. Input Dataset

The environmental forcing data used to drive the drift simulations were sourced from the Copernicus Marine Service (CMEMS), one of the six thematic services of the European Union’s Copernicus Earth Observation Programme. CMEMS was selected as the primary data provider for its ability to deliver high-resolution, quality-controlled oceanographic and atmospheric products that are continuously updated through the assimilation of satellite and in situ observations, ensuring the accuracy and temporal consistency essential for reliable drift modelling. Moreover, the free and open access policy of the service, combined with the broad spatial coverage of its datasets, makes it a well-established reference in the operational oceanography and SAR communities.

Specifically, the extracted data include: Stokes Drift obtained from the Mediterranean Sea Forecasting System; wind data and sea surface current provided by the global ocean datasets.

The Stokes drift was extracted from the Mediterranean Sea Waves Analysis and Forecast product, which provides hourly wave parameters at a horizontal resolution of 1/24° (approximately 4.2 km) [16,17]. The dataset is built on the WAM 4.6.2 wave model and structured as a nested sequence of two computational grids to ensure that swell propagating from the North Atlantic is correctly channelled through the Strait of Gibraltar into the Mediterranean basin. The dataset includes integrated parameters such as significant wave height, period, direction, and the zonal and meridional components of Stokes drift.

Wind data were instead drawn from the global Level-4 hourly product at a spatial resolution of 0.125° (approximately 14 km), which combines scatterometer observations from Metop-B and Metop-C ASCAT with collocated ECMWF operational model fields to correct for persistent biases in the hourly numerical output, ultimately providing stress-equivalent wind variables together with their divergence and curl [10,18].

Finally, sea surface currents were obtained from the Global Ocean Physics Analysis and Forecast product at a horizontal resolution of 1/24° (approximately 4.2 km) [19,20]. The model component is based on the NEMO platform forced at the surface by ECMWF ERA5 reanalysis fields, with observations assimilated through a reduced-order Kalman filter that jointly ingests along-track altimetric sea level anomaly data, satellite sea surface temperature, sea ice concentration, and in situ vertical profiles of temperature and salinity.

2.4. Simulated Scenarios

Under ideal conditions, characterised by a perfectly known flow field and a fully comprehensive numerical model, it would be possible to reconstruct with precision the trajectories of the bodies recovered along the shore and to identify the exact location of the shipwreck. In practice, however, both the input datasets and the numerical models are subject to inherent limitations that prevent such a deterministic reconstruction. These limitations stem from the difficulty of accurately representing the full range of natural processes that affect the drift, including unresolved small-scale dynamics, atmospheric and oceanic variability, and uncertainties in the initial conditions. Therefore, under ideal conditions, it would be sufficient to release a single particle for each recovery point to reconstruct the trajectory accurately. In practical applications, however, this approach is not feasible, and a different strategy must be adopted to properly account for the uncertainties in the input data and the inherent variability of the environmental conditions.

The simulations were carried out by initialising particle releases at the recovery sites and propagating them backwards in time within the numerical model. For each recovered body, both the location and the date of recovery were used as initial conditions for the simulations. The Leeway module was configured to simulate the drift of Person-In-Water (PIW) – deceased/incapacitated objects (object type 6). To account for the uncertainties associated with both the recovery location and the exact time at which the bodies reached the shore, each source was represented by an ensemble of 1000 drifters, uniformly distributed within a circular area of 1000 m radius. For a given source, the ensemble was released with temporal offsets relative to midnight of the recovery day. Specifically, five different release times were considered: 0 h, 12 h, 24 h, 48 h, and 72 h. The release time of the drifters was estimated as midnight of the recovery day minus the assigned temporal offset. Each combination of recovery site and release time constituted an independent simulation scenario, resulting in a total of 25 scenarios (5 bodies × 5 release times), each comprising 1000 drifter trajectories. All simulations were integrated backwards from the date of recovery back to 1 February 2024, resulting in simulation durations of approximately 45 days for ID 1 (recovered on 17 March 2024), 72 days for ID 2 (recovered on 13 April 2024), 73 days for ID 3 (recovered on 14 April 2024), 77 days for ID 4 (recovered on 18 April 2024), and 150 days for ID 5 (recovered on 30 June 2024).

2.5. Trajectory-Based Probability Analysis and Shipwreck Reconstruction

A systematic analysis of the resulting drifter trajectories was performed in order to identify the most probable origin areas of the recovered bodies. The study area was discretised into a regular grid of approximately equal-sized cells, over which the trajectory data were spatially aggregated.

For a given reference time $t_r$, each grid cell $c_i$ was analysed by counting the number of simulated drifter trajectories intersecting the cell within a prescribed temporal tolerance $\Delta t$ centred on $t_r$. Only trajectories satisfying both spatial and temporal consistency criteria were retained in the analysis. Let $N_i\left(t_r\right)$ denote the number of drifters passing through cell $c_i$ at time $t_r$, and let $S_i\left(t_r\right)$ be the number of distinct recovery locations contributing to those trajectories. The likelihood of cell $c_i$ being the origin of the recovered remains can be expressed as:

$$P_d(c_i,t_r)=D_t·D_d={\displaystyle \frac{N_i\left(t_r\right)}{{max}_j{N}_j\left(t_r\right)}}·{\displaystyle \frac{S_i\left(t_r\right)}{N_{\mathrm{source}}}}$$

(5)

where $N_{\mathrm{source}}$ is the total number of recovery locations. The first term ($D_t$) normalises the local drifter density by the maximum value observed across all cells, providing a relative measure of trajectory concentration. The second term ($D_d$) accounts for the number of independent recovery sources whose trajectories pass through the cell, thereby penalising areas reached by drifters from only a subset of the recovery locations. This dual formulation ensures that high-likelihood areas reflect not only a dense convergence of trajectories but also a consistent spatial overlap across multiple independent sources, thus reducing the influence of localised or spurious pathways that may arise from a single recovery event.

To provide a quantitative measure of the degree of hydrodynamic connectivity between the study domain and the recovery locations at each reference time $t_r$, a scalar index, hereinafter referred to as the Hydrodynamic Connectivity Index (HCI), was introduced. The index is defined as the product of two complementary components: the mean drift probability over the active cells, representing the intensity of the connectivity, and the fractional area of the domain exhibiting significant connectivity, representing its spatial extent. Formally:

$${\mathrm{HCI}}_c\left(t_r\right)={\displaystyle \frac{{\sum }_i{P}_d(c_i,t_r)·\mathbf{1}[P_d(c_i,t_r)>\theta ]}{{\sum }_i\mathbf{1}[P_d(c_i,t_r)>\theta ]}}·{\displaystyle \frac{N_{\theta }\left(t_r\right)}{N_{\mathrm{tot}}}}$$

(6)

where $\theta$ is a prescribed probability threshold, $\mathbf{1}[·]$ is the indicator function, $N_{\theta }\left(t_r\right)$ is the number of grid cells whose drift probability exceeds $\theta$ at time $t_r$, and $N_{\mathrm{tot}}$ is the total number of marine grid cells in the study domain. The first term represents the mean value of $P_d$ over the subset of cells exhibiting significant connectivity (i.e., those with $P_d>\theta$), providing a measure of the intensity of the hydrodynamic pathway. The second term quantifies the fractional area of the domain involved in the connectivity, penalising configurations in which high probability values are confined to a small number of isolated cells. The HCI thus attains high values only when elevated drift probabilities are simultaneously widespread and intense across the domain, making it a robust indicator of the establishment of a coherent large-scale transport pathway towards the recovery locations. For the present analysis, the threshold was set to $\theta =0.05$. The index was subsequently normalised by its maximum value over the analysis period to facilitate temporal comparison:

$$HCI\left(t_r\right)={\displaystyle \frac{{\mathrm{HCI}}_c\left(t_r\right)}{{max}_{t_r}\mathrm{HCI}\left(t_r\right)}}$$

(7)

The temporal evolution of $HCI$ provides a compact scalar summary of the spatiotemporal development of the hydrodynamic connectivity identified by the backward drift simulations.

The resulting probability maps were subsequently cross-referenced with documented shipwreck incidents reported in the central Mediterranean between January and February 2024, retrieved from publicly available institutional and journalistic sources. Once this correspondence was established, a dedicated analysis of the wave climate was carried out over the area encompassing the most probable transit corridor of the vessel. The analysis aimed to assess whether severe sea state conditions consistent with the capsizing of a structurally inadequate vessel were present along the reconstructed route during the estimated time window of the incident, in order either to support or to rule out the hypotheses regarding the vessels’ routes.

3. Results

For the post-processing and analysis of the simulation results, a specific spatial domain was defined within the geographical boundaries of $7.{50}^{\circ }\mathrm{E}$ to ${15}^{\circ }\mathrm{E}$ longitude and ${35}^{\circ }\mathrm{N}$ to ${45}^{\circ }\mathrm{N}$ latitude. To compute spatial statistics and aggregate trajectory behaviours, this domain was discretised into a regular grid composed of cells with a spatial extent of $30\mathrm{km}$.

The month of February and part of March were designated as the temporal reference period for the analysis. Within this framework, the estimation of the drift probability ($P_d$) was conducted using a computational time step ($\Delta t$) of $2.5\mathrm{day}$, with reference dates established every five days.

Figure 3 and Figure 4 show the spatial and temporal evolution of the drift source density ($D_d$) and drift trajectory density ($D_t$) over the study area, evaluated at five-day intervals from 1 February to 12 March 2024.

Figure 5 shows the spatial and temporal evolution of the drift probability $P_d$ over the study area. Each panel represents a snapshot of the hydrodynamic connectivity between a given grid cell and the observed recovery locations at the corresponding reference time $t_r$. Specifically, $P_d(c_i,t_r)$ is interpreted here as the probability that an object located in cell $c_i$ at time $t_r$ would have experienced, over the subsequent temporal evolution of the flow field, the oceanographic and meteorological conditions necessary to be transported towards the recovery sites.

In the earliest panels (1–11 February 2024), the probability field is largely uniform and low across the entire domain, indicating that no spatially coherent pathway connecting the study area to the recovery locations had yet been established. Starting from approximately 13 February, a region of elevated $P_d$ begins to emerge at latitudes around 37°–38°N and longitudes between 8° and 11°E, corresponding to the area of the eastern Strait of Sicily and the Tunisian coastal waters.

To contextualise the results of the backward drift simulations within the broader framework of migration events occurring in the study period, a systematic review of publicly available sources was conducted to identify shipwreck incidents reported in the Central Mediterranean between January and February 2024 (see Table 2).

The review identified four documented incidents in the central Mediterranean route during the period of interest (Table 2). On 11 January 2024, a boat carrying 37 Tunisian nationals – predominantly young people aged between 13 and 35 from the village of El Hencha in the Sfax Governorate – departed from Sfax at approximately 14:30 local time; all contact with the vessel was lost by 22:00 the same day, and the boat has never been found [21]. On 5–6 February 2024, a vessel carrying at least 17 Tunisian nationals, including a five-year-old child, departed from the port of Bizerte and subsequently disappeared in the central Mediterranean with no survivors found [22]. On 5–6 February 2024, a rubber dinghy carrying 55 migrants of African nationalities departed from Al-Zawiya, Libya, at approximately 23:00 local time and capsized approximately six hours later after taking on water north of the coastal city of Zuwara, leaving 53 people dead or missing, including two babies; only two survivors were rescued [23]. On 8 February 2024, a metal boat carrying 42 Sudanese migrants capsized nine nautical miles off the port of Chebba, resulting in 13 confirmed deaths and 27 missing persons, with only two survivors rescued by the Tunisian coast guard [24].

Among the incidents identified, only the Bizerte shipwreck of 5–6 February 2024 is spatially and temporally compatible with the probability maps derived from the backward drift simulations. The port of Bizerte is located at the northernmost tip of Tunisia (approximately 37.27°N, 9.87°E), and vessels departing from this location bound for Italy would have followed a north-northeastward trajectory across the central Mediterranean, consistent with the high-probability region identified in the $P_d$ maps between mid-February and early March 2024. The departure date of 5-6 February 2024 falls within the temporal window estimated by the simulations, and Bizerte’s geographic position lies within the area of elevated hydrodynamic connectivity identified in the analysis. The remaining three incidents, all originating from the Sfax–Chebba coastal stretch (approximately 34°–35°N, 11°E), are located approximately 400 km southeast of the high-probability region and are therefore not consistent with the reconstructed drift pathways. Based on this analysis, the Bizerte shipwreck of 5–6 February 2024 is identified as the most probable event associated with the recovery of the five bodies examined in the present study.

To further investigate the probable causes of the shipwreck, wave data were extracted at selected grid points within the area that may have been transited by the vessel during its crossing, based on the high-probability regions identified by the $P_d$ analysis and the known departure point of Bizerte. Specifically, time series of significant wave height and associated spectral parameters were retrieved from the CMEMS dataset at grid points corresponding to the most probable transit corridor of the embarkation, to assess the sea state conditions prevailing along the reconstructed route at the time of the estimated incident. The location of the grid points is shown in Figure 2. The rationale behind this analysis lies in the well-documented relationship between severe sea states and shipwreck events: overloaded and structurally inadequate vessels, such as those typically employed in irregular migration crossings, are particularly vulnerable to capsizing under storm wave conditions, where significant wave heights and steep wave periods can rapidly exceed the stability thresholds of such embarkations. The occurrence of intense storm events along the estimated transit corridor during the temporal window identified by the backward drift simulations was therefore examined as an independent line of evidence to corroborate the spatiotemporal reconstruction and to provide a physically consistent explanation for the vessel’s sinking.

Figure 6 shows the time evolution of the significant wave height extracted at the five locations (see Figure 2) and the value of HCI estimated for each reference date.

The wave data reveal that the study area was affected by several distinct storm events throughout the period, with significant wave height ($H_{m0}$) values exceeding 3–4 m at multiple grid points. A first significant storm event is clearly identifiable around 11–12 February 2024, when $H_{m0}$ reached a peak of approximately 4.9 m at point P3, with all five grid points recording values above 2 m simultaneously. A second, more prolonged period of elevated wave heights is observed between late February and early March 2024, with $H_{m0}$ values consistently exceeding 2–3 m across all points, reaching a peak of approximately 5.5 m around 4–5 March 2024. The red numbers annotated along the top axis denote the reference dates of the corresponding $P_d$ probability maps shown in Figure 5, facilitating direct temporal comparison between the wave climate and the evolution of the hydrodynamic connectivity.

The wave analysis alone, however, is insufficient to uniquely identify the location or timing of the shipwreck. The occurrence of multiple severe storm events throughout the entire period and across the full spatial extent of the study area implies that wave conditions alone cannot discriminate between competing hypotheses regarding the departure point or the site of the incident. This limitation underscores the importance of integrating the wave data with the results of the backward drift simulations.

Specifically, when all available lines of evidence are considered jointly, the spatial distribution of $P_d$ and its temporal evolution as quantified by HCI, the geographic and temporal compatibility with the documented shipwreck incidents, as well as the severe sea state conditions prevailing along the estimated transit corridor — the most probable scenario is that the five bodies recovered in the Aeolian Islands area between March and June 2024 were among the passengers of the vessel that departed from Bizerte on 5-6 February 2024. However, the progressive increase in HCI observed from mid-February onwards, peaking around 4–5 March 2024, is not consistent with the hypothesis that the sinking occurred during or shortly after the first intense storm event of 11–12 February 2024, and most likely the vessel encountered damage after the departure, thus causing its shipwreck.

4. Discussion and Conclusions

The present study proposes and evaluates a methodology for reconstructing the most probable origin and trajectory of human remains recovered at sea, through the combination of large-scale numerical model outputs — including meteorological, oceanographic, and wave forecasting datasets — with Lagrangian drift simulations. The approach was applied to a real case study involving five bodies recovered in the Aeolian Islands area between March and June 2024, which represents an ideal test bed for assessing the reliability and limitations of the proposed framework.

A key aspect of the methodological design adopted in the present study is that the analysis was conducted under the assumption that no prior information regarding the departure point of the vessel was available. This deliberate choice was made in order to test the autonomous predictive capability of the methodology, independently of the forensic and investigative findings that subsequently identified Bizerte (Tunisia) as the most probable point of departure. Under this blind hypothesis, backward Lagrangian simulations were carried out using environmental forcing data provided by the Copernicus Marine Service in combination with the OpenDrift Leeway model [11]. To account for the inherent uncertainties associated with recovery location, time of death, and environmental forcing, a multi-scenario simulation framework was adopted, releasing ensembles of 1000 drifters per scenario from each recovery site across five temporal offsets (0, 12, 24, 48, and 72 hours), resulting in a total of 25 independent simulation scenarios.

The results of the backward drift simulations were synthesised through the drift probability metric $P_d$, which quantifies the hydrodynamic connectivity between each grid cell of the study domain and the observed recovery locations at each reference time. The temporal evolution of $P_d$ was further condensed into the scalar Hydrodynamic Connectivity Index (HCI), which combines the intensity and spatial extent of the high-probability regions into a single diagnostic quantity. The analysis revealed a progressive emergence and intensification of a coherent high-probability region in the area between the Tunisian coastal waters and the Strait of Sicily, centred at approximately 37°–38°N and 8°–10°E, starting from mid-February 2024 and reaching its maximum extent between late February and early March 2024. This region is geographically consistent with the Central Mediterranean migration corridor connecting the northern Tunisian coast to the Aeolian Islands area.

The cross-referencing of the $P_d$ maps with documented shipwreck incidents reported in the central Mediterranean during January–February 2024 identified the Bizerte shipwreck of 5-6 February 2024 as the only event spatially and temporally compatible with the simulation results. This finding is in full agreement with the conclusions independently reached by the forensic investigations conducted on the recovered bodies, thereby validating the proposed methodology. The complementary analysis of the wave climate along the estimated transit corridor led to the observation of severe storm events — with significant wave heights exceeding 4–5 m — starting from 11-12 February, not consistent with the incident, excluding that the capsizing of the vessel could be due to severe storm events.

The proposed approach presents several strengths that make it suitable for operational application in forensic SAR contexts. It is entirely based on freely accessible, high-resolution operational datasets; it is applicable when multiple bodies are recovered from different locations and at different times, as the multi-source formulation of $P_d$ explicitly exploits the spatial consistency across independent recovery events to reduce ambiguity; and it provides a probabilistic, quantitative framework that can be readily integrated with forensic, pathological, and investigative evidence.

Nevertheless, the methodology is subject to a number of limitations that should be acknowledged. The accuracy of the reconstructed trajectories is inherently constrained by the spatial and temporal resolution of the large-scale numerical models employed, which may not adequately resolve mesoscale oceanographic features or localised atmospheric forcing. Furthermore, the Leeway model relies on empirically derived coefficients that may not fully capture the hydrodynamic behaviour of the specific type of vessel involved. The temporal window of the simulations, while informed by forensic estimates of time of immersion based on Lepas spp. colonisation, remains subject to biological variability. Finally, the approach is most effective when applied to cases involving multiple recovery events, as the spatial consistency criterion embedded in the $P_d$ formulation requires at least two independent sources to provide meaningful discrimination between competing origin hypotheses.

In conclusion, the results of the present study demonstrate that the integration of backward Lagrangian drift simulations with large-scale operational oceanographic and atmospheric datasets constitutes a viable and reliable tool for the forensic reconstruction of shipwreck events in the central Mediterranean. The methodology successfully identified the Bizerte departure as the most probable origin of the five recovered individuals, in agreement with independent forensic findings, and provided a physically consistent reconstruction of the conditions that likely led to the sinking of the vessel.