Assessing lagoon plume dispersion under variable metocean conditions

Highlights

• Plume dispersion analyzed under different metocean scenarios.
• Particle exposure maps reveal accumulation areas along the lagoon’s boundary.
• Periodic release simulations show that plume dispersion sensitivity to release timing is advection-driven.
• Mixing properties highlight variations in transport mechanisms across scenarios.

1 Introduction

Coastal lagoons are among the most dynamic and ecologically significant environments. Despite occupying a small fraction of the global coastal ocean, these semi-enclosed water bodies sustain high biodiversity and provide critical habitat and nursery functions for many species. Their restricted connection to the open sea, however, makes them particularly vulnerable to anthropogenic pressures associated with densely populated and economically active coastlines, including urban development, agriculture, fisheries, and tourism (Neumann et al., 2015). As a consequence, lagoons and other semi-enclosed basins increasingly experience water-quality deterioration, eutrophication, habitat degradation, and resource overexploitation (Raicevich et al., 2018; Rodrigues et al., 2021; Soria et al., 2022).

From a physical standpoint, lagoon circulation is strongly controlled by the limited water exchange with the adjacent shelf and by local forcing mechanisms. Wind stress, tides, and the geometry and management of inlets and hydraulic structures can modulate transport pathways, thereby influencing where and how rapidly nutrients, contaminants, and suspended matter accumulate (Rynne et al., 2016; Simonetti and Cappietti, 2022). Because these controls operate on short spatial scales and are highly sensitive to small morphological or operational changes, predicting dispersion and exposure patterns is often essential for environmental management and risk mitigation.

To support monitoring and protection of lagoon ecosystems, two complementary strategies are typically adopted: (i) observations, increasingly integrating remote sensing and automated analytics, and (ii) numerical modelling of hydrodynamics and transport. Traditional field campaigns rely on in situ measurements of key physico-chemical parameters (e.g., temperature, dissolved oxygen, chemical oxygen demand, nutrients) to quantify water quality and diagnose eutrophication (Scanes et al., 2007; Sigua et al., 2020). More recently, satellite observations (e.g., Sentinel-2 and Landsat) combined with machine learning techniques have enabled spatially extensive and frequent retrieval of turbidity and chlorophyll-related proxies, supporting near-real-time assessments and, in some cases, operational monitoring frameworks (Li et al., 2020; Erena et al., 2019; Caballero et al., 2022; Mohan et al., 2025). In parallel, modelling studies commonly investigate lagoon circulation and lagoon–sea exchange to estimate renewal metrics, accumulation zones and to reproduce the spatio-temporal distribution of water-quality variables (Soria et al., 2021; Torres Bejarano et al., 2020; Ingrassia et al., 2024; Zannella et al., 2025).

This study investigates the dispersion of floating passive tracers using a set of pre-computed weekly scenarios designed to represent the metocean variability of the study domain. For each scenario, we analyse Lagrangian trajectories through statistical descriptors commonly used to characterise transport and mixing in geophysical flows (LaCasce, 2008). In addition, we quantify the sensitivity of plume evolution to release timing by comparing an instantaneous release with a temporally distributed release, separating the role of advective pathway variability from plume spreading.

The Narta Lagoon (southern Albania) provides an ideal testbed for this analysis: it is shallow, ecologically valuable, and characterised by limited and intermittently managed exchange with the Adriatic Sea. Previous study (Scovenna et al., 2025) revealed how weak is the exchange of water between the lagoon and the outer sea. For this reason, here the lagoon is treated as closed basin. Particle releases are prescribed at a location, identified as the most contaminated area of the lagoon through previous sampling analysis (Bakaj et al., 2022). The objective is therefore to simulate a real-world discharge and its behavior within a closed system, using virtual particles as a first approximation. Results show that different weekly circulation regimes can lead to markedly different plume behaviours, ranging from persistent confinement near the source to enhanced lagoon-wide spreading under more mixing-favourable conditions. The paper is organised as follows: Section 2 describes the study area; Section 3 presents the modelling framework and analysis methods; results are reported in Section 4 and discussed in Section 5, along with limitations and future improvements; conclusions are summarised in Section 6.

2 Study area

The Narta Lagoon lies on the southern Albanian coast near the city of Vlorë (Figure 1, Panels A,–B). The lagoon is recognized for its ecological and socio-economic relevance, as it provides habitat for migratory birds and European eels and supports local aquaculture and related activities that employ a substantial workforce in the area (Elmazi and Bazini, 2011; Kane et al., 2015; Mladenov et al., 2017; Topi et al., 2013).

In recent years, the system has been exposed to multiple anthropogenic stressors. Potential sources of pressure include waste disposal facilities serving Vlorë and nearby towns, industrial activities such as salt production and leather tanning, and illegal fishing practices (Kane et al., 2015; Çako et al., 2014; Beqaj, 2023; Vaso et al., 1996; Topi et al., 2013). Several of these activities occur close to the shoreline, and incomplete wastewater treatment raises concerns about pollutant inputs via surface runoff and groundwater pathways.

Water exchange with the Adriatic Sea occurs through two inlets with contrasting morphology and management (bottom panel in Figure 1). The southern channel (19.395647°E, 40.514683°N; WGS84) is ∼ 200 m long, with a minimum width of ∼ 10 m, and includes manually operated sluice gates at its narrowest section. The northern channel (19.379781°E, 40.545369°N) is longer (∼ 800 m) and slightly wider (> 10 m), but it is prone to seasonal closure due to sedimentation driven by littoral transport; local reports indicate that it is manually reopened by fishermen roughly twice per year (Topi et al., 2013). Overall, inlet operation is largely informal, with no fixed schedule for opening or maintenance, which can promote stagnation and ecological stress. Moreover, even under idealized conditions with fully open inlets, modelling results indicate weak lagoon–sea exchange, with the external water level exerting a negligible influence on the lagoon circulation except in the immediate vicinity of the inlets (Scovenna et al., 2025). Based on this evidence, the lagoon is treated here as a closed basin, to focus on the internal transport mechanisms shaping the distribution of virtual particles.

3 Methodology

To drive the Particle Tracking Model (PTM), six weekly surface velocity fields are extracted from a high-resolution 3D hydrodynamic model run with Delft3D 4 (Lesser et al., 2004). Given the assumption that the lagoon is a closed basin, the model is only forced with wind, namely, 6 weekly scenarios representative of a 40-year historical record.

Using the surface layer of the velocity fields, the PTM simulates the dispersion of passive tracers by numerically integrating particle trajectories under advection and diffusion:

x t + Δ t = x t + ∫ t t + Δ t u x , t d t + ∇ K Δ t + R 2 K * Δ t (1)

where x is the particle position, u ( x , t ) is the velocity at that position (obtained by linear interpolation of the surrounding grid velocities), and K is the horizontal diffusion coefficient. Because K varies spatially (see Equation 2), a naive random walk based on the stochastic term alone would violate the well-mixed condition, causing spurious accumulation of particles in low-diffusivity regions. To prevent this, a gradient term ∇ K Δ t is added explicitly to Equation 1, following Visser (1997), and the stochastic term is evaluated using K *, the diffusivity at the midpoint position reached after applying half of the gradient drift, which improves the accuracy of the correction. R = ( R x , R y ) consists of two independent random numbers, one per spatial component, each drawn from a Gaussian distribution with zero mean and unit variance, N ( 0,1 ), to ensure isotropic diffusion. Δ t is the integration time step and integration is performed using a fourth-order Runge–Kutta scheme, while, following the Smagorinsky approach, the diffusion coefficient K is parameterized in terms of the local shear of the flow:

K = C s Δ x Δ y ∂ u ∂ x 2 + ∂ v ∂ y 2 + 1 2 ∂ u ∂ y + ∂ v ∂ x 2 , (2)

where Δ x Δ y is the grid-cell area and C s is the Smagorinsky coefficient. Although the Smagorinsky parameterization was originally formulated as a closure for subgrid-scale momentum, it is commonly repurposed to estimate horizontal eddy diffusivity in Lagrangian particle-tracking applications (van Sebille et al., 2018), implicitly assuming a turbulent Schmidt number of unity ( S c t = ν t / K h = 1 ). This simplification is becoming popular in particle-tracking studies of pollutant and microplastic transport (Mishra et al., 2024), where direct estimates of tracer diffusivity are typically unavailable and the velocity-shear-based closure offers a practical, spatially-varying alternative to a constant diffusivity. While S c t may deviate from unity in stratified or strongly anisotropic flows, we adopt this standard assumption following common practice in the field, and the magnitude of K obtained is compared against the Lagrangian eddy diffusivity inferred independently from particle trajectory statistics (Section 4.3) as a consistency check.

Post-processing first quantifies exposure, defined here as the time-averaged particle occupancy within each model cell, normalized over a specific temporal window. This window, denoted as T plateau, corresponds to the period during which the number of active particles released reaches and maintains a plateau. Let n i j ( t n ) denote the number of particles located inside cell ( i , j ) at the hourly time t n and let N plateau bet the total number of time steps within the plateau window T plateau. The time-averaged occupancy in cell ( i , j ) during this window is defined as:

n ̄ i j = 1 N plateau ∑ t n T plateau n i j t n , (3)

Single-particle statistics are then computed to characterise dispersion scales and to compare scenarios, namely, the Lagrangian velocity autocorrelation (LaCasce, 2008):

R τ = u t ⋅ u t + τ ‖ u t ‖ 2 , (4)

where ⟨ ⋅ ⟩ denotes an ensemble average over particles, and the absolute dispersion, i.e., the mean-squared displacement from the release position (LaCasce, 2008). Let x i ( t ) be the position vector of particle i at time t and N the number of particles in the considered ensemble. The absolute dispersion is defined as

A 2 t = 1 N ∑ i = 1 N x i t − x i 0 2 , (5)

which combines both the drift of the cloud and its spreading. Here, along with presenting A 2 as a time series, we also map it in 2D by binning particles according to their release positions, to visualise spatial patterns in displacement. The absolute dispersion is further decomposed into the coherent drift, quantified through

A c 2 t = x ̄ t − x ̄ 0 2 , (6)

where x ̄ ( t ) = 1 N ∑ i = 1 N x i ( t ) represents the ensemble centroid, and the spreading around the drifting centroid, quantified by the relative dispersion,

R c 2 t = 1 N ∑ i = 1 N x i ′ t 2 , (7)

where x i ′ ( t ) = x i ( t ) − x ̄ ( t ).

Finally, the anisotropy is computed using the eigenvalues of the 2 × 2 spatial covariance tensor for each particle cohort. The components of the covariance tensor are given by:

Σ t = σ x x t σ x y t σ x y t σ y y t , (8)

with the individual terms calculated as the ensemble average over the active trajectories ( N ) within a cohort:

σ x x t = 1 N ∑ i = 1 N x i ′ t 2 , σ y y t = 1 N ∑ i = 1 N y i ′ t 2 , σ x y t = 1 N ∑ i = 1 N x i ′ t y i ′ t . (9)

By performing an eigen decomposition of Σ ( t ), we obtain the eigenvalues λ 1 ( t ) and λ 2 ( t ), representing the variances along the major and minor principal axes of the particle cloud, respectively. The axis-independent anisotropy is then assessed using:

A N t = λ 1 t − λ 2 t λ 1 t + λ 2 t . (10)

The orientation of the principal axis of stretching relative to the original coordinate system is determined by the eigenvector [ v x , v y ] T associated with the major eigenvalue λ 1:

θ t = atan 2 v y , v x , (11)

where θ ( t ) is the principal axis angle in radians. Values of A N → 1 indicate highly directional, line-like stretching along the angle θ, while A N = 0 indicates perfectly isotropic, circular spreading.

4 Results

4.1 Clustering and hydrodynamics

The principal characteristics of the hydrodynamic model are summarized in Table 1. The computational domain is discretized using a regular horizontal grid with a spatial resolution of 30 m and 7 terrain-following σ-layers in the vertical dimension. To ensure numerical stability and maintain a Courant number below 5, the integration time step was set to 9 s. Additionally, a 12-h spin-up period was implemented to allow the flow field to adjust to the imposed initial and boundary conditions. Vertical turbulent mixing was parameterized using the standard k − ϵ closure scheme. The model is forced by six spatially uniform, weekly wind scenarios extracted from a 40-year hindcast produced by the MeteOcean research group at the University of Genoa (Lira-Loarca et al., 2022; Cassola et al., 2015; Mentaschi et al., 2013). This hindcast provides hourly data on a grid with a spatial resolution of approximately 10 km covering the entire Mediterranean Sea. The representative scenarios were selected by applying the Maximum Dissimilarity Algorithm (MDA) (Camus et al., 2011) to the wind hindcast at the node closest to the study area. This node is located inside the lagoon at the position marked by the white ‘X’ in Figure 1. Further details regarding the scenario selection methodology are available in Scovenna et al. (2025), and the corresponding wind roses are presented in Figure 2. The resulting weekly averaged surface velocity fields are shown in Figure 3; Scenarios 01 and 02 (top panels) exhibit overall weak velocities, whereas Scenario 06 is the most intense, featuring westward-directed currents that align coherently with the dominant wind forcing.

TABLE 1. Hydrodynamic model characteristics.

4.2 Dispersion simulation

For the Lagrangian simulations, the integration time step is set to Δ t = 60 s, with particle positions recorded every 30 min. The Smagorinsky coefficient is parameterized as C s = 0.1, yielding a horizontal diffusion coefficient K that ensures turbulent spreading without introducing excessive artificial diffusion (see Table 2 in Section 4.3 for a quantitative assessment of K). Sensitivity analyses indicated that higher values of C s introduce numerical artifacts near boundaries, where inflated diffusive steps can cause near-shore particles to artificially cross the land-water interface, thereby neutralizing the physical effects of the slip conditions described hereafter.

TABLE 2. Comparison of the prescribed Smagorinsky diffusivity ( K ̄ ) and the Lagrangian eddy diffusivity ( κ iso ) derived from Taylor’s (1921) theory, for each scenario.

To prevent artificial beaching, a slip condition is applied at the boundaries as a velocity of 0.01 m/s pointing inside the lagoon: whenever a particle approaches the shore, it is gently moved away from it, along the local normal direction.

The initial release positions of the particles is selected based on environmental sampling analysis from previous studies (Bakaj et al., 2022), which identified the area adjacent to a wastewater discharge channel as the most contaminated zone. The geometry of the release area is configured to mimic a localized plume, represented by the gray shaded region in Figure 5.

A sensitivity analysis is conducted to determine the optimum number of particles required for statistical convergence in the simulation. The same simulation is iteratively run releasing each time a different amount of particles. For each run, the final spatial distribution of particles is mapped onto the hydrodynamic grid and normalized by the total number of particles released. These normalized concentration fields is then compared against a reference simulation by calculating the Root Mean Squared Error (RMSE) across the domain. The convergence threshold is identified where an increase in particle number yields diminishing returns in RMSE reduction. The reference simulation utilizes 46,494 particles released at a 5.3-m spacing. As shown in Figure 4, the black arrow indicates the selected optimal particle count. The figure legend details the particle count and corresponding initial spacing for each run; for example, the chosen point represents a scenario with n p = 9184 particles released at a 13-m spacing.

A total of n p virtual particles are therefore uniformly released at the sewage outfall, under two configurations. In the first configuration, all particles are released at the initial time step. In the second configuration, ensembles of n p particles are released every 6 h during the first 92 h of the simulation to assess the sensitivity of plume dispersion to release timing. To ensure each cohort remains in the simulation for the same amount of time, each particle is removed 72 h after its release. Because historical information on discharge magnitude and composition is unavailable, release characteristics are prescribed heuristically.

In the first set of simulations, particles are released only at the initial time step (instantaneous release). Figure 5 shows the relative trajectories across all scenarios. The grey cloud in the bottom-right corner of the lagoon indicates the release area, while the red dots mark the final particle positions. Intermediate positions are coloured according to simulation time, illustrating the progressive dispersion of the plume. Marked differences emerge across scenarios. In scenario 01, particles remain clustered near the release area; in scenarios 02, 03 and 04 particles occupy the southern half of the lagoon, whereas in scenario 05, particles are mainly trapped by boundary-following currents. Finally, in scenario 06 particles are moved toward the western part of the lagoon and kept close to the shore for the rest of the simulation.

Figures 6, 7 show the exposure maps for the second configuration (continuous release) across all simulated scenarios. The exposure is computed by applying Equation 3 with a temporal window of T plateau = 24 h. To assess the sensitivity of the system to boundary dynamics, the occupancy is computed under both no-slip (left panels) and slip (right panels) conditions; the former considers particles to be immediately beached upon entering a land cell.

A logarithmic color scale is implemented to enhance the visibility of low-occupancy areas, while red boxes denote cells with the highest occupancy (top 0.5%) to highlight primary accumulation zones. Notably, the discrepancy between the slip and no-slip conditions is minimal, demonstrating that the identification of these accumulation areas is remarkably robust against the choice of boundary treatment. This indicates that particle behavior and retention close to the shore are fundamentally driven by the local hydrodynamics themselves, rather than being an artifact of the numerical slip condition.

Consequently, both configurations yield identical localized features: the high-occupancy clusters consistently locate within the southern half of the lagoon, accumulating either along the eastern shore (S01, S02, S05) or the western shore (S03, S04, S06), where, even within the analyzed closed-system configuration, particles exhibit a strong tendency to accumulate near the northern inlet (see Figure 1). If properly managed, this inlet could therefore offer a viable pathway for particle flushing and export from the system.

4.3 Mixing properties

Starting from the instantaneous-release configuration, Figure 8 shows the Lagrangian velocity autocorrelation (Equation 4) for the u (blue) and v (orange) components; lighter lines correspond to individual trajectories, while darker lines indicate the ensemble mean. A marked difference emerges between the scenario 01 and the others. In S01, the autocorrelation remains persistent and exhibits a clear oscillatory pattern, particularly for the v component, indicating the presence of coherent and/or periodically forced motions. Notably, the thin lines closely follow the average signal, suggesting particles tend to move as a whole, rather than as a collection of independent trajectories. By contrast, in the others scenarios the autocorrelation decays rapidly toward zero, suggesting a much shorter Lagrangian memory time scale and a more rapidly decorrelating flow, while, except for S05, the individual particle signals are much more dispersed around the mean: particles quickly lose memory of their initial velocity and start behave independently.

The Lagrangian integral time scale T L is computed as the integral of the velocity autocorrelation function from 0 to its first zero crossing, and is reported as a vertical dashed line on the absolute dispersion plots (Figure 9). It marks the expected transition between the ballistic regime (A 2 ∝ t 2 for t < T L) and the diffusive regime (A 2 ∝ t for t > T L), whose slopes are shown as black dashed reference lines. Thin orange lines represent A 2 (Equation 5) for individual particles, while the bold orange line is the ensemble mean. The ballistic regime is evident in all scenarios; however, beyond T L the absolute dispersion does not transition to the expected diffusive growth, instead plateauing (e.g., in S01) or exhibiting sharp oscillatory increases and decreases (e.g., in S02). The latter behavior is attributed to oscillatory motion driven by the unsteady, wind-driven hydrodynamics. In both cases, the diffusive regime expected for stationary, homogeneous turbulence is masked by coherent advective transport.

The Taylor-derived eddy diffusivity, κ iso = σ u 2 T L, where σ u 2 is the variance of the Lagrangian velocity fluctuations (Taylor, 1921), is computed for each scenario and compared against the time- and space-averaged Smagorinsky diffusivity(Equation 2). Results are summarized in Table 2. In all scenarios, κ iso / K ̄ ≫ 1, confirming that total particle dispersion is governed by resolved-scale advection rather than by the parameterized sub-grid diffusion, consistent with the absence of a clear diffusive regime in the A 2 curves.

Figure 10 displays the spatial distribution of the absolute dispersion at t = 168 h. Here, A 2 is normalized between 0 and 1 to emphasize behavioral differences within the plume rather than absolute magnitude. In most scenarios, the normalized A 2 is approximately uniform across the plume, indicating that dispersion is largely insensitive to release position. Scenario 02 is a notable exception: particles released in the upper half of the plume tend to disperse farther from their release position, while those released in the lower half remain comparatively closer.

Mixing descriptors are also applied to the continuous-release experiment. Results relative to S01 and S03 are presented, as representative of the two main distinct behaviors, while results relative to the other scenarios can be found in the Supplementary Materials. Figure 11 compares S01 (left column) and S03 (right column), where each colour denotes a different release cohort (from light blue to violet the cohorts are released at subsequent times). The analysis for different scenarios yields different y-axis values, as the focus is not on the magnitude differences between the descriptors across scenarios, but rather on the behavior of different cohorts within the same scenario. The first row shows the absolute dispersion A 2 ( t ). In S01, cohorts exhibit broadly similar behaviour, indicating limited sensitivity to release time. In S03, by contrast, cohorts released earlier (light blue colours) differ substantially from those released later (violet colours), highlighting a pronounced dependence on release timing. The A 2 ( t ) curves exhibit highly variable initial values, revealing that as the simulation progresses, the current velocity increases and, consequently, cohorts released later are immediately advected far from the releasing position. This behavior is also evident in the light blue curves, which exhibit a mild initial slope followed by a sharp, steep growth after approximately 12 h.

To disentangle whether these differences originate from coherent transport or from plume spreading, the second and third rows report the squared displacement of the center of mass, A c 2 ( t ) (Equation 6), and the dispersion relative to the center of mass, R c 2 ( t ) (Equation 7), respectively. Consistent with A 2 ( t ), S01 shows only minor inter-cohort variability in both drift and spreading. In S03, the cohort-to-cohort differences are primarily expressed in the drift term A c 2 ( t ), whereas the spreading term R c 2 ( t ) remains comparatively similar across cohorts. This indicates that the sensitivity to release time in S03 is mainly controlled by changes in advective pathways rather than by changes in the intensity of spreading, which results in particles of the same cohort being transported coherently.

Figure 12 displays the anisotropy metric (top panel, computed applying Equations 8–10) alongside the orientation of the principal axis of deformation (bottom panel, Equation 11), which has been mapped symmetrically to enhance visualization. These metrics reveal markedly distinct hydrodynamic behaviors between the two scenarios. In S01, there is no dominant or persistent deformation direction; the anisotropy index continuously oscillates across a wide range ( 0 ≤ A N ≤ 0.8 ), indicative of intermittent isotropic spreading. Conversely, in S03, all cohorts exhibit highly directional stretching, maintaining a near-constant A N ∼ 0.9 throughout most of the week. The corresponding rose plots corroborate these findings, demonstrating a highly constrained alignment for S03, where the principal axis of deformation is sharply clustered around 150 °.

5 Discussion

To investigate the dispersion of floating virtual particles released near a sewage outfall, we simulated the resulting trajectories under six weekly metocean scenarios, using two release configurations: a single release at the initial time, and a continuous release in which a new cohort was introduced every 6 h and tracked for 72 h. Trajectories from the continuous-release configuration were analyzed through exposure maps, while single-particle dispersion metrics (autocorrelation and absolute dispersion) were computed for both configurations.

Exposure maps were computed with and without a slip condition at the lagoon boundaries; the two formulations yielded no significant differences, with particles consistently accumulating along the lagoon margins. Notably, even though the lagoon is treated as a closed system, particles tend to accumulate near the north inlet in several scenarios, underscoring how important properly managing this inlet would be in allowing particles to exit the lagoon.

For the single-release configuration, both the autocorrelation and the absolute dispersion show markedly different behavior between low- and high-energy scenarios. Under low-energy conditions—S01 in particular—particles remain tightly clustered and move largely as a single cloud, staying close to the release point. Under high-energy conditions, stronger mixing disperses the particles more uniformly across the lagoon. The absolute dispersion time series exhibits a ballistic regime up to the Lagrangian decorrelation time, after which no diffusive regime emerges; this is consistent with the flow being advection-dominated rather than diffusion-dominated beyond that timescale. The spatial distribution of absolute dispersion shows coherent behavior within the plume across all scenarios except S02, where particles released in the upper half of the plume disperse markedly farther from their release position than the rest of the cohort. Comparing the continuous-release results for S01 and S03 shows that absolute dispersion is largely insensitive to release timing in S01. In S03, by contrast, the sensitivity to release timing is driven primarily by the advective component of the flow: varying the release time does not change how far particles spread around their center of mass, but does change which advective pathway the cloud follows.

Several limitations should be acknowledged. In the absence of historical release data, the number of particles and the release area were chosen heuristically, so the present setup should be regarded as a test case rather than a site-calibrated configuration; using historical data to constrain these parameters would improve the realism of the results. Likewise, the use of virtual, passive particles—rather than tracers with the decay, settling, or buoyancy properties of specific contaminants—limits how directly these results translate into pollutant-dispersion risk assessments. The lack of in situ data precluded any calibration of the Lagrangian simulations; future work should calibrate the model against drifter observations. Finally, given the results of the exposure analysis, showing particles accumulating close to the north inlet, it would be worthy to run the same Lagrangian simulation with the open-inlets hydrodynamics, in order to quantify the percentage of particles leaving the lagoon.

6 Conclusion

This study examines the dispersion of virtual particles released near a sewage outfall in the Narta Lagoon under six weekly metocean scenarios, using two particle-release configurations. The results highlight how hydrodynamic conditions and release timing jointly shape the distribution and accumulation of particles within the plume. In several scenarios, dispersion proved highly sensitive to release timing, with the advective component of the flow playing the dominant role in this variability. The analysis further demonstrates the value of exposure and dispersion descriptors for characterizing the spatial dynamics of floating particles in semi-enclosed systems: by simulating a range of hydrodynamic conditions and release configurations, we identified the areas most susceptible to particle accumulation and gained insight into the lagoon’s mixing properties.

Taken together, these results contribute to the growing body of knowledge on lagoon dynamics and pollutant dispersion, and offer a framework that can support the management of this sensitive ecosystem. More broadly, the scenario-based approach is flexible and computationally efficient, and can be extended to other systems where limited exchange with the open sea poses similar challenges for water-quality management. Further refinement of the model, together with the inclusion of pollutant-specific properties, will improve its predictive capacity and support better-informed decision-making for the protection of these fragile ecosystems.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenve.2026.1882143/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source

Footnotes

Edited by: Ciro Fernando Bustillo LeCompte, Toronto Metropolitan University, Canada

Reviewed by: Saeed Hariri, Institut Français de Recherche pour l'Exploitation de la Mer, France

Davide Pasquali, University of L'Aquila, Italy

Martin F. Soto-Jimenez, National Autonomous University of Mexico, Mexico