ChemicalDrift 1.0: an open-source Lagrangian chemical-fate and transport model for organic aquatic pollutants

A new model for transport and fate of chemicals in the aquatic environment is presented. The tool, named ChemicalDrift, is integrated into the open-source Lagrangian framework OpenDrift and is hereby presented for organic compounds. The supported chemical processes include the degradation, the volatilization, and the partitioning between the different phases that a target chemical can be associated with in the aquatic environment, e.g. dissolved, bound to suspended particles, or deposited to the seabed sediments. The dependencies of the chemical processes on changes in temperature, salinity, and particle concentration are formulated and implemented. The chemical-fate modelling is combined with wide support for hydrodynamics by the integration within the Lagrangian framework which provides e.g. advection by ocean currents, diffusion, wind-induced turbulent mixing, and Stokes drift generated by waves. A flexible interface compatible with a wide range of available metocean data is made accessible by the integration, making the tool easily adaptable to different spatio-temporal scales and fit for modelling of complex coastal regions. Further inherent capabilities of the Lagrangian approach include the seamless tracking and separation of multiple sources, e.g. pollutants emitted from ships or from rivers or water treatment plants. Specific interfaces to a dataset produced by a model of emissions from shipping and to an unstructured-grid oceanographic model of the Adriatic Sea are provided. The model includes a database of chemical parameters for a set of poly-aromatic hydrocarbons and a database of emission factors for different chemicals found in discharged waters from sulfur emission abatement systems in marine vessels. A post-processing tool for generating mean concentrations of a target chemical, over customizable spatio-temporal grids, is provided. Model development and simulation results demonstrating the functionalities of the model are presented, while tuning of parameters, validation, and reporting of numerical results are planned as future activities. The ChemicalDrift model flexibility, functionalities, and potential are demonstrated through a selection of examples, introducing the model as a freely available and open-source tool for chemical fate and transport that can be applied to assess the risks of contamination by organic pollutants in the aquatic environment.

1 Introduction

The negative effects of chemical pollution on the environment and human health have long been established (Naidu et al., 2021), but only more recently has the need for integrated strategies and solutions for the simultaneous management of multiple stressors been widely recognized (Pirotta et al., 2022). In fact, one of the main challenges in assessing the impact of anthropogenic stressors on the environment, and especially on aquatic ecosystems, is to estimate the exposure to these contaminants while also keeping track of their origin, in particular when multiple sources of chemical pollution are present.

Lagrangian models (van Sebille et al., 2018) offer the possibility of investigating the transport of chemicals without losing information on their origin, thus facilitating the assessment of the contribution each source has to the overall exposure of the target ecosystem. Furthermore, since their spatial resolution is not inherently limited to the resolution of Eulerian grid boxes and no grid generation is required (Azarpira et al., 2021), Lagrangian models are suitable for studying areas with complex geometries such as coasts, lagoons and archipelagos. The capabilities of these models to predict the transport of chemicals and particles emitted into the atmosphere by a variety of sources (e.g. industrial installations, Hirdman et al., 2010; Lee et al., 2014, nuclear accidents, Becker et al., 2007; Draxler et al., 2015, and volcanic eruptions, Prata et al., 2007; D'Amours et al., 2010) have been widely studied (Onink et al., 2021). On the other hand, the potential of Lagrangian models applied to the fate and transport of organic chemicals in the water compartment is not fully investigated.

It is recognized that coastal environments are of particular concern due to the presence of multiple stressors (e.g. wastewater discharge, soil and sediment contamination, agricultural and municipal run-off and combustion of fossil fuels by civil and industrial activities) (Danovaro and Boero, 2019; Szymczycha et al., 2019; Calgaro et al., 2019, 2021), posing a risk to the ecosystem and human health due to water and air contamination. An additional specific stressor is represented by emissions of SO_x, NO_x and other airborne contaminants (e.g. particulate matter and polyaromatic hydrocarbons) from shipping activity, which has been shown to significantly degrade air quality (Zhang et al., 2021), especially in areas characterized by high traffic such as ports and shipping lanes. For these reasons, the International Maritime Organization (IMO), starting from 1 January 2020, has introduced new global regulations limiting the sulfur content of fossil fuels to a maximum of 0.50  % weight / weight ($w/w$), unless they are fitted with a suitable exhaust-gas-cleaning system (“scrubbers”) (Vedachalam et al., 2022). Furthermore, the more stringent limit of 0.10  % $w/w$ sulfur has been applied to especially sensitive areas (i.e. emission control areas, ECAs) like the Baltic Sea, the North Sea and the US coasts (Chang et al., 2018). Plans to include other areas in the ECA list, like the Mediterranean Sea, are under discussion (Testa, 2020).

The scrubbing process is implemented by leading the exhaust gas through a fine-water mist where SO_x and NO_x are dissolved together with particulate matter, metals and organic compounds (Lunde Hermansson et al., 2021), and thus large volumes of toxic (Turner et al., 2017) and highly acidified effluents (Turner et al., 2018) (where several pollutants, e.g. Cu, Pb, Hg, V, Zn, As, Ni and polyaromatic hydrocarbons (PAHs), have been detected) are produced and discharged into the sea. Moreover, due to the increase in shipping traffic predicted for the next decades (Sardain et al., 2019), it is crucial to develop methods and tools to investigate the impact of scrubber use on the marine ecosystem in relation to other sources of chemical pollution.

The present work is carried out under the scope of the Horizon 2020 EMERGE project, which aims to develop an integrated modelling framework to assess the combined impact of shipping emission-control options on the aquatic and atmospheric environment. Utilizing the open-source Lagrangian framework OpenDrift (Dagestad et al., 2018), a new model for chemical transport and fate of aquatic pollutants, named ChemicalDrift, has been developed and is hereby presented and demonstrated.

This article is organized as follows: the relevant chemical processes implemented in ChemicalDrift are presented and formulated in Sect. 2; the set of inputs utilized in this work, including metocean forcing data, chemical parameters for the studied pollutants, and emission source data, are presented in Sect. 3, and a selection of demonstrative examples are presented in Sect. 4. Planned future developments and conclusions are summarized in Sect. 5.

2 Chemical-fate modelling

The new module ChemicalDrift is presented in this section, supported by a detailed description of the modelled chemical processes. The module is coded in Python and integrated into the open-source Lagrangian framework OpenDrift (Dagestad et al., 2018), where the fundamental physical processes of advection and diffusion as well as atmospheric and oceanographic forcing were previously implemented. The OpenDrift framework already contains a module for oil pollution, OpenOil (Dagestad et al., 2018), which has been validated and used for simulations of oil spills (Röhrs et al., 2018; Brekke et al., 2021; Hole et al., 2021), i.e. mixtures of petroleum products modelled as oil droplets or oil/water emulsions in the water column. The module presented here has a different target, being intended for simulations of dispersion of specific single chemicals, which might be hazardous contaminants even at very low concentrations in the water column or in the sediment layer, and includes a new implementation of the relevant chemical processes.

Only the chemistry of non-ionizable organic compounds is described in this work, since the chemistry of metals has already been implemented in the Radionuclides OpenDrift module (Simonsen et al., 2019a), which is integrated into ChemicalDrift.

The reactions implemented in ChemicalDrift include a partitioning scheme between the different phases that a target chemical can be associated with within the aquatic environment, the degradation of organic chemicals in the water column and in the sediments and the volatilization of dissolved chemicals from the water to the atmosphere. The dependencies of the each reaction on temperature and salinity changes are described and implemented based on the scientific literature. Sedimentation and resuspension are also presented in this section since these physical processes are strongly related to the chemical-phase partitioning.

2.1 Dynamic partitioning between compartments

Chemicals in the aquatic environment can be associated with different components of this medium (i.e. dissolved in the water or bound to dissolved organic carbon (DOC), suspended particle matter (SPM) or sediments) and are thereby exposed to different physical and chemical processes. For example, dissolved chemicals will be transported by advection due to water currents and by turbulent diffusion processes, such as wind-induced vertical mixing, while chemicals adsorbed to solid particles will also be affected by gravity and might sink towards the sea floor. Furthermore, the different distributions between dissolved and bound chemicals will also influence its availability to degradation processes like hydrolysis, biodegradation and photolysis. In computational models it is therefore crucial to process each of these components in distinct compartments.

In ChemicalDrift, each Lagrangian element represents a given mass of a target chemical, and the partitioning between the media components is implemented by assigning each element to a single corresponding model compartment. The partitioning is carried out dynamically, as opposed to a simpler steady-state approach, and transfer rates between compartments are calculated and applied at each time step. The dynamic approach is fit for the changing environmental conditions each Lagrangian element is exposed to (i.e. water temperature, salinity, SPM and DOC concentration), as it is transported in the media.

The dynamic partitioning scheme used in ChemicalDrift was first introduced by Periáñez and Elliott (2002) and implemented in the OpenDrift Radionuclides module (Simonsen et al., 2019a, b), where it was applied to metals. In this work, several adaptations were made to expand its applicability to organic chemicals and consider the influence of temperature and salinity changes on the transfer rates.

The main phases considered in the algorithm (Fig. 1) include (i) dissolved in the water fraction, (ii) adsorbed to the SPM fraction, and (iii) the sediment fraction. The sediment fraction includes both the chemical adsorbed to solid particles that have settled to the sea floor and the chemical dissolved in the interstitial water present in the space between the settled particles. Optional model compartments can also be activated. The chemical can also be modelled as adsorbed to DOC and form other aggregates that will be subjected to sinking (i.e. the DOC phase). Adsorption to solid particles can be treated as a two-step process where the diffusion of the chemical within the particle is also considered (i.e. a slowly reversible fraction, SR) for both SPM and sediments. Finally, an optional compartment for buried sediments is implemented, representing a sink for the target chemical that is progressively buried below layers of overlying sediments and is therefore not available for exchanges with the water column.

Figure 1 The model compartments implemented in ChemicalDrift and the corresponding exchange processes and transfer rates adapted from Simonsen et al. (2019a). The default (dissolved, adsorbed to SPM and sediments) compartments are indicated with solid contours. Optional compartments are depicted with dashed contours.

2.1.1 Adsorption and desorption rates

Transfer rates are calculated by adapting the equations described by Simonsen et al. (2019a) and are applied using the approach explained in Sect. 2.1.2 for estimating the adsorption and desorption kinetic constants, k_ads (L kg⁻¹ h⁻¹) and k_des (h⁻¹).

Transfer rates from the dissolved phase to the DOC, SPM and sediment compartments (i.e. k₁₂, k₁₃ and k₁₄) and the corresponding desorption rates (i.e. k₂₁, k₃₁ and k₄₁) are given by

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}& k_{\mathrm{12}}=k_{\mathrm{ads}}{C}_{\mathrm{DOC}},\\ & k_{\mathrm{13}}=k_{\mathrm{ads}}{C}_{\mathrm{SPM}},\\ & k_{\mathrm{14}}=k_{\mathrm{ads}}{\mathit{\rho }}_{\mathrm{sed}}(\mathrm{1}-p_{\mathrm{sed}})L\mathit{\varphi }\mathit{\delta }/H,\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}& k_{\mathrm{21}}=k_{\mathrm{des}}=k_{\mathrm{ads}}/K_{d,\mathrm{DOC}},\\ & k_{\mathrm{31}}=k_{\mathrm{des}}=k_{\mathrm{ads}}/K_{d,\mathrm{SPM}},\\ & k_{\mathrm{41}}=\mathit{\varphi }{k}_{\mathrm{des}}=\mathit{\varphi }{k}_{\mathrm{ads}}/K_{d,\mathrm{SPM}}.\end{array}\end{array}$$

Here the indices refer to the numbering illustrated in Fig. 1: C_DOC is the concentration of DOC in the water column (kg L⁻¹), C_SPM is the concentration of SPM in the water column (kg L⁻¹), ρ_sed is the density of sediment solids (kg m⁻³), p_sed is the porosity of the sediments, L (m) is thickness of the sediment layer which interacts with the water column, H (m) gives the thickness of the water layer above the sea bottom that interacts with the sediment, ϕ is the fraction of effective sorbents considering that a fraction of the sediment particle surface may be unavailable for sorption–desorption interactions due to neighbouring sediment particles and δ is the logical variable which ensures that only the chemicals in the seabed interaction layer are able to adsorb to sediments and is set to 1 when the distance to the seabed is less than H and 0 elsewhere. Default values or typical ranges are summarized in Table 1. The concentrations C_DOC and C_SPM are provided as modelled or experimentally derived data fields. If only 2-D surface data are available, the vertical profiles are estimated assuming a constant concentration in the mixed layer (equal to the surface values) and exponentially decreasing with depth below the mixed layer, with user-definable constants, in accordance with the observations of biogeochemical Argo float measurements reported by Galí et al. (2022).

Fantke (2018); Honti et al. (2016); Kettelarij et al. (2019)Paterson and Mackay (1995); Kettelarij et al. (2019)Periáñez (2012)Kobayashi et al. (2007); Periáñez and Elliott (2002)Periáñez (2008); Periáñez et al. (2013)

Table 1 Parameters for the calculation of sediment-layer adsorption and desorption rates.

The stochastic method for calculating the probability of phase change between model compartments given the time step and the transfer rates is described by Periáñez and Elliott (2002).

2.1.2 General estimation of k_ads and k_des for organic chemicals

The partitioning of organic chemicals to organic carbon has been reported as involving rapid adsorption to particle surfaces, followed by slow movement into, and out of, organic matter and porous aggregates (Karickhoff and Morris, 1985), thus determining the time needed for the attainment of equilibrium (Park and Clough, 2014).

The solid–water partition coefficient K_d (L kg⁻¹) is the ratio between the adsorption and desorption kinetic constants:

$$\begin{array}{}\text{(3)}& K_d=k_{\mathrm{ads}}/k_{\mathrm{des}}.\end{array}$$

Given this context, Karickhoff and Morris (Karickhoff and Morris, 1985; Site, 2001) found the reciprocal of k_des to be a linear function of K_d over 3 orders of magnitude (r²=0.87), following $\mathrm{1}/k_{\mathrm{des}}\approx \mathrm{0.03}{K}_d$. Using this approximation and the definition of K_d, the adsorption and desorption kinetic constants can be estimated as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& k_{\mathrm{ads}}\approx \mathrm{33.3}\mathrm{L}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}},\\ & k_{\mathrm{des}}=k_{\mathrm{ads}}/K_d\approx \mathrm{33.3}/K_d{\mathrm{h}}^{-\mathrm{1}}.\end{array}\end{array}$$

The value of K_d for a given chemical in the aquatic environment can be subjected to significant changes due to e.g. temperature, salinity, SPM and DOC concentration variations and pH; therefore, the use of experimental local measurements of K_d is recommended. However, if such data are not available, the model will estimate the K_d as described below.

For organic chemicals it is assumed that the organic carbon fraction of solid particles (i.e. suspended particle matter and sediments) is almost entirely responsible for their sorbing capacity, and therefore K_d can be expressed as the product of a water–organic carbon partition coefficient (i.e. K_OC – L kg${}_{\mathrm{OC}}^{-\mathrm{1}}$) and the organic carbon content of the dry matter (f_OC – g_OC g⁻¹).

$$\begin{array}{}\text{(5)}& K_d=K_{\mathrm{OC}}{f}_{\mathrm{OC}}\end{array}$$

Typical values for f_OC are in the range 0.01–0.10 g_OC g⁻¹ (Seiter et al., 2004).

Relationships between K_OC and the octanol–water partition coefficient K_OW (L kg⁻¹) are available from previous studies. In the case of non-ionizing organic chemicals, the following formulations reported by Park and Clough (2014) are used to estimate K_OC for SPM and DOC, respectively.

$$\begin{array}{}\text{(6)}& {\displaystyle }{K}_{\mathrm{OC},\mathrm{DOC}}\approx \mathrm{2.88}(K_{\mathrm{OW}}{)}^{\mathrm{0.67}}\text{(7)}& {\displaystyle }{K}_{\mathrm{OC},\mathrm{SPM}}\approx \mathrm{2.62}(K_{\mathrm{OW}}{)}^{\mathrm{0.82}}\end{array}$$

2.1.3 Effect of temperature and salinity on the solid–water partition coefficient

Experimental values of K_OC and K_d are usually obtained at standard or ambient temperature and for fresh water, but temperature and salinity variations have been shown to have a strong impact on the adsorption of chemicals to SPM, DOC, and sediments. The following equation and correction factors are applied to obtain to obtain the solid–water partition coefficient K_d(T,S) for a desired value of temperature T (K) and salinity S (PSU) based on a known reference value K_d(T_ref,0):

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}{K}_d(T,S)& =K_d(T_{\mathrm{ref}},\mathrm{0})\mathrm{TempCorr}(T,T_{\mathrm{ref}},\mathrm{\Delta }H)\\ & {\mathrm{SalinityCorr}}_{K_d}\left(S\right),\end{array}\text{(9)}& {\displaystyle }{\mathrm{SalinityCorr}}_{K_d}\left(S\right)={\mathrm{10}}^{k_{\mathrm{S}}{C}_{\mathrm{salt}}\left(S\right)},\text{(10)}& {\displaystyle }\mathrm{TempCorr}(T,T_{\mathrm{ref}},\mathrm{\Delta }H)=e^{-\frac{\mathrm{\Delta }H}{R}\left(\frac{\mathrm{1}}{T}-\frac{\mathrm{1}}{T_{\mathrm{ref}}}\right)},\end{array}$$

where k_S is the Setschenow (salting-out) constant in seawater of the non-ionizable organic chemical with respect to the specific salt considered (L mol⁻¹), and C_salt is the molar concentration of the salts in solution (mol L⁻¹). Salt concentration is calculated starting from salinity considering the molecular weight of seawater salt, MolWt_salt=68.35 g mol⁻¹ (Schwarzenbach et al., 2016), and the water density, ρ.

$$\begin{array}{}\text{(11)}& C_{\mathrm{salt}}=\left({\displaystyle \frac{S}{{\mathrm{MolWt}}_{\mathrm{salt}}}}\right)\mathit{\rho }\end{array}$$

The correction factor for temperature is based on the Arrhenius equation where ΔH (J mol⁻¹) is the enthalpy of adsorption for the given chemical and R is the universal gas constant (R=8.3145 J mol⁻¹ K⁻¹). The model utilizes distinct enthalpies for DOC and SPM, indicated as ΔH_KOC,DOC and ΔH_KOC,SPM, respectively.

2.1.4 Sedimentation and resuspension

Particle sedimentation and resuspension dynamics in ChemicalDrift are partly based on methods previously implemented in other OpenDrift modules. The vertical motion of Lagrangian elements representing dissolved chemicals or chemicals adsorbed to DOC is calculated as for a passive tracer by adding vertical currents (typically small compared to the other terms) to the vertical mixing associated with wind-induced or background turbulent diffusion.

Lagrangian elements associated with SPM are also affected by gravity and will sink towards the sea floor, and hence a third component has to be added to calculate the vertical motion. This term is referred to as the terminal velocity and is calculated from Stokes' law (Stokes, 1851), assuming the particles to be spherical with known density and diameter (Sundby, 1983). Empirical relationships are utilized to update the water density (Fofonoff and Millard Jr, 1983) and viscosity (Ådlandsvik, 2000; Myksvoll et al., 2013) based on T and S as the elements sink.

When sinking particles reach the sea floor, the chemical elements are transferred to the sediment-layer compartments and thus can be subjected to either sediment burial or resuspension. Resuspension occurs when the horizontal current at the bottom is greater than a given threshold v_crit. Resuspended elements are lifted from the previous sea-floor depth z_min (depth below sea level, negative) to a depth $z_{\min }+z_{\mathrm{res}}+\mathcal{N}\left(z_{\mathrm{0},\mathrm{res},\mathrm{unc}}\right)$, where z_res is the resuspension depth, z_res,unc is the resuspension depth uncertainty and $\mathcal{N}(\mathrm{0},z_{\mathrm{res},\mathrm{unc}})$ is a randomly generated number from a normal distribution with zero mean and standard deviation z_res,unc (Simonsen et al., 2019a).

2.2 Degradation

Organic chemicals in the aquatic environment can be degraded due to various processes, such as biodegradation, photodegradation and hydrolysis. Modelling each of these reactions separately requires a large amount of information on the environmental behaviour of the target chemical and on the characteristics of the selected study area. Since these data are typically difficult to obtain with sufficient accuracy for a wide range of chemicals and case-study regions, a simpler approach is implemented in the current version of ChemicalDrift: degradation is implemented considering distinct overall reaction rates for the water column (subscript Wtot) and the sediments (subscript Stot) and modelled as a first-order kinetic decay of the mass associated with each Lagrangian element; in detail, overall degradation rate constants k or the corresponding half-lives $t_{\mathrm{1}/\mathrm{2}}$ can be obtained experimentally without separating the effects of the three sub-processes. The mass m_i(t) of each Lagrangian chemical element is calculated at each time step applying

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& m_i(t+\mathrm{\Delta }t)=\\ & \left\{\begin{array}{ll}{m}_i\left(t\right)e^{-k_{\text{W\hspace{0.17em}tot}}\mathrm{\Delta }t},& \text{if }i\in (\mathrm{dissolved},\mathrm{DOC}),\\ m_i\left(t\right)e^{-k_{\text{S\hspace{0.17em}tot}}\mathrm{\Delta }t},& \text{if }i\in \left(\mathrm{sediments}\right),\\ m_i\left(t\right),& \text{otherwise },\end{array}\right.\end{array}\end{array}$$

where k_W tot is the overall (total) degradation constant for the chemical either dissolved in the water column or adsorbed to dissolved organic carbon/matter, k_S tot is the overall degradation constant for chemicals in the sediment layer and Δt is the time-step interval. As reported by, among others, Kettelarij et al. (2019), it is assumed that degradation of chemicals adsorbed to suspended particles can be neglected. Degradation rate constants can be specified directly or calculated from the corresponding half-life constants $t_{\mathrm{1}/\mathrm{2}\mathrm{W}\mathrm{tot}}$, $t_{\mathrm{1}/\mathrm{2}\mathrm{S}\mathrm{tot}}$ applying $k=-\log \left(\mathrm{0.5}\right)/t_{\mathrm{1}/\mathrm{2}}$. Furthermore, the effect of temperature is considered by applying to each degradation rate constant a correction factor calculated by applying the Arrhenius equation as in Eq. (10) using the reference temperature (i.e. T_{ref W tot} and T_{ref S tot}) and enthalpy (i.e. ΔH_W tot and ΔH_S tot) associated with each process. Salinity has also been demonstrated to have an impact on the microbial populations responsible for the biodegradation of organic chemicals, but since the modelling of these effects would require the integration of ChemicalDrift with a biogeochemical model, it is considered outside the scope of this work and is neglected (Park and Clough, 2014).

2.3 Volatilization

Volatilization is modelled using the “stagnant boundary theory”, or two-film model, in which the target chemical must diffuse across both a stagnant water layer and a stagnant air layer to volatilize out of the water column (Schwarzenbach et al., 2016; Parnis and Mackay, 2020; Park and Clough, 2014; Kettelarij et al., 2019).

According to the literature (Schwarzenbach et al., 2016), the volatilization rate (kg m⁻³ s⁻¹) can be calculated based on the properties of the two-film interface and the concentrations of the chemical:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathrm{volatilization}& ={\displaystyle \frac{\partial }{\partial t}}{C}_{\mathrm{water}}=\\ & -({\mathrm{MTC}}_{\mathrm{vol}}/H_{\mathrm{MLD}})(C_{\mathrm{sat}}-C_{\mathrm{water}}).\end{array}\end{array}$$

Here, C_water (kg m⁻³) is the concentration of the chemical in the water phase, C_sat (kg m⁻³) is the saturation concentration of the pollutant in equilibrium with the gas phase, H_MLD (m) is the thickness of the mixed layer if the water column is stratified or the maximum depth otherwise and MTC_vol (m s⁻¹) is the overall volatilization mass transfer coefficient from water to air.

If pollutants can be assumed to have a negligible concentration in air (i.e. C_sat≈0), then the equation simplifies to a first-order problem (i.e. $\frac{\partial }{\partial t}{C}_{\mathrm{water}}=k_{\mathrm{vol}}{C}_{\mathrm{water}}$) that is solved by $C_{\mathrm{water}}\left(t_{\mathrm{2}}\right)=C_{\mathrm{water}}\left(t_{\mathrm{1}}\right)e^{-k_{\mathrm{vol}}(t_{\mathrm{2}}-t_{\mathrm{1}})}$, where the volatilization rate constant (s⁻¹) can be calculated from Eq. (14).

$$\begin{array}{}\text{(14)}& k_{\mathrm{vol}}={\mathrm{MTC}}_{\mathrm{vol}}/H_{\mathrm{MLD}}\end{array}$$

This assumption is often considered acceptable and is widely applied for chemical-fate modelling (Park and Clough, 2014; Deltares, 2023), especially in the case of chemicals with very low volatility or that can be found mostly bound to suspended particles in the air, such as the five- and six-ring PAHs, but they can be problematic for more volatile chemicals (e.g. naphthalene) that are found in significant concentrations in the gas phase. In those cases, a coupled atmosphere–ocean model is required for accurate calculation of exchanges across the interface, but ChemicalDrift can offer a compromise between accuracy of results and model complexity.

Table 2 Chemical parameters utilized for modelling the sorption, degradation, and volatilization of three organic compounds.

In the Lagrangian framework utilized by the proposed model, concentrations are not calculated at each time step, and hence the equivalent effect is obtained by applying an exponential decay to the mass m_i(t) of all dissolved elements i within the mixed layer, as expressed by $m_i(t+\mathrm{\Delta }t)=m_i\left(t\right)e^{-k_{\mathrm{vol}}\mathrm{\Delta }t}$.

As reported above, k_vol can be estimated starting from MTC_vol, whose reciprocal can be interpreted as the total resistance to the mass transfer through the two-film interface and therefore can be expressed as the sum of the resistances of each layer (Eq. 15):

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\mathrm{1}}{{\mathrm{MTC}}_{\mathrm{vol}}}}={\displaystyle \frac{\mathrm{1}}{{\mathrm{MTC}}_{\mathrm{w}}}}+{\displaystyle \frac{\mathrm{1}}{{\mathrm{MTC}}_{\mathrm{a}}{H}^{\mathrm{cc}}}}.\end{array}$$

Here, MTC_a and MTC_w (m s⁻¹) are the air-side and water-side mass transfer coefficients, respectively, and H^cc is the dimensionless Henry law constant.

The diffusion rates of a chemical in these stagnant boundary layers can be related to the known diffusion rates of reference substances such as oxygen and water vapour (Schwarzenbach et al., 2016). Based on the boundary-layer theory (Deacon, 1977), the mass transfer coefficient for a target chemical can be obtained from these by applying a correction factor that depends on the ratios of the Schmidt numbers of the target and reference chemicals. The following relations proposed by McGillis et al. (2001) and (Johnson, 2010) and summarized by Schwarzenbach et al. (2016) are used in this work.

$$\begin{array}{}\text{(16)}& {\displaystyle }\begin{array}{rl}& {\mathrm{MTC}}_{\mathrm{a}}=\mathrm{0.01}{\left({\displaystyle \frac{{\mathrm{MolWt}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}{\mathrm{MolWt}}}\right)}^{\mathrm{1}/\mathrm{3}}{\mathrm{MTC}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{Oa}}\\ & {\mathrm{MTC}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{Oa}}=\mathrm{0.1}+\\ & {\displaystyle \frac{U_{\mathrm{10}}(\mathrm{6.1}+\mathrm{0.63}{U}_{\mathrm{10}}{)}^{\mathrm{1}/\mathrm{2}}}{\begin{array}{c}\mathrm{13.3}(S{c}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{Oa}}{)}^{\mathrm{1}/\mathrm{2}}+(\mathrm{6.1}\times {\mathrm{10}}^{-\mathrm{4}}+\mathrm{6.3}\times {\mathrm{10}}^{-\mathrm{5}}{U}_{\mathrm{10}}{)}^{\mathrm{-}\mathrm{1}/\mathrm{2}}\\ \mathrm{-}\mathrm{5}+\mathrm{1.25}\ln \left(S{c}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{Oa}}\right)\end{array}}}\end{array}\text{(17)}& {\displaystyle }\begin{array}{rl}{\mathrm{MTC}}_{\mathrm{w}}& =\mathrm{0.01}(\mathrm{9}\times {\mathrm{10}}^{\mathrm{-}\mathrm{4}}+\mathrm{7.2}\times {\mathrm{10}}^{\mathrm{-}\mathrm{6}}{U}_{\mathrm{10}}^{\mathrm{3}})\\ & {\left({\displaystyle \frac{{\mathrm{MolWt}}_{{\mathrm{CO}}_{\mathrm{2}}}}{\mathrm{MolWt}}}\right)}^{\mathrm{1}/\mathrm{4}}\end{array}\end{array}$$

Here, MolWt (g mol⁻¹) is the molecular weight of the target chemical, ${\mathrm{MolWt}}_{{\mathrm{CO}}_{\mathrm{2}}}$ is the molecular weight of CO₂ (44 g mol⁻¹), ${\mathrm{MolWt}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ is the molecular weight of water (18 g mol⁻¹), U₁₀ is the wind speed 10 m above the water surface (m s⁻¹) and ${\mathit{Sc}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{Oa}}$ is the Schmidt number of water vapour in air (0.62). Then, Eqs. (16) and (17) are inserted into Eq. (15), where the dependencies on temperature and salinity are accounted for when the dimensionless Henry law volatility constant H^cc is calculated (Park and Clough, 2014) as

$$\begin{array}{}\text{(18)}& H^{\mathrm{cc}}={\displaystyle \frac{H^{\mathrm{pc}}\left(\mathrm{1}+\mathrm{0.01143}S\right)}{RT}}.\end{array}$$

Here, R is the universal gas constant ($\mathrm{8.206}\times {\mathrm{10}}^{-\mathrm{5}}$ atm m³ mol⁻¹ K⁻¹), T is the temperature (K) and S is the salinity (PSU). The last term at the numerator is a dimensionless linear-correction factor for the effect of salinity, accounting for a factor of 1.4 for S=35 PSU compared to fresh water (Park and Clough, 2014; Schwarzenbach et al., 1993). H^pc is Henry's volatility constant (atm m³ mol⁻¹), which is either specified by the user or estimated by

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& H^{\mathrm{pc}}=\\ & {\displaystyle \frac{P^{\mathrm{sat}}\left(T_{\mathrm{ref},\mathrm{vp}}\right)\mathrm{TempCorr}(T,T_{\mathrm{ref},\mathrm{vp}},\mathrm{\Delta }{H}_{\mathrm{vol}})}{C_W^{\mathrm{sat}}\left(T_{\mathrm{ref},\mathrm{sol}}\right)\mathrm{TempCorr}(T,T_{\mathrm{ref},\mathrm{sol}},\mathrm{\Delta }{H}_{\mathrm{sol}})}}\mathrm{MolWt},\end{array}\end{array}$$

where P^sat (atm) is the vapour pressure of the target chemical at the reference temperature, C^sat (g m⁻³) is the solubility of the target chemical at the reference temperature, and $\mathrm{TempCorr}(\cdot ,\cdot ,\cdot )$ are dimensionless temperature correction factors calculated with Eq. (10) for volatilization and solubilization reactions with the respective enthalpy changes ΔH_vol and ΔH_sol at the target temperature T.

Table 3 CMEMS products and data variables.

3 Emissions, chemical parameters and metocean forcing

3.1 Ship Traffic Emission Assessment Model (STEAM)

Emission fields from STEAM by Jalkanen et al. (2021) are used in this work. Modelled emissions from open-loop scrubbers in 2019 are utilized in the following examples. The method named seed_from_STEAM(…) is implemented in the ChemicalDrift class to allow seamless seeding of Lagrangian elements directly from the STEAM DataArray. Emission data are provided as a volume (m³) of discharged scrubber water at each spatio-temporal grid point. The conversion from volume to mass of chemicals is implemented within the class method, applying average emission factors (µg L⁻¹) and confidence intervals reported in Lunde Hermansson et al. (2021) for open_loop and closed_loop scrubbers for a set of heavy metals and PAHs, selected with the arguments scrubber_type and chemical_compound, respectively. The total mass of chemicals for each STEAM data point is then distributed on a number of Lagrangian elements with a given mass and over a circular area defined by the arguments mass_element_ug (µg) and radius (m).

3.2 Chemical parameters for PAHs

A database of parameters for a set of organic compounds and heavy metals has been compiled under the scope of the EMERGE project, based on an extensive review of data available in the literature. Mean values from the database for a selection of PAHs are integrated into the ChemicalDrift class method init_chemical_compound, which allows one to configure the whole set of parameters by simply selecting the target chemical. Alternatively, each parameter can be specified singularly with the set_config method. Both procedures are shown in the example in Listing 1, while the data for the PAHs utilized in this work are reported in Table 2.

Listing 1 Chemical simulation set-up and configuration.

3.3 Metocean forcing from CMEMS and the System of HydrodYnamic Finite Element Modules (SHYFEM)

Meteorological and oceanographic forcing is obtained from different sources. Several CMEMS products are utilized in this work, providing water temperature, salinity, current horizontal velocities, mixed-layer depth, ocean surface winds, and SPM concentration, as summarized in Table 3.

Listing 2 Seeding elements and running simulation.

High-resolution 3-D currents, water temperature and salinity forcing over the northern Adriatic Sea including the Venice Lagoon are provided by the application of the unstructured hydrodynamic SHYFEM model (Bellafiore et al., 2018). SHYFEM model integrations are performed on a spatial domain covering the whole of the Adriatic Sea and its main coastal lagoons (Marano-Grado, Venice, and the Po River delta). To adequately resolve the river–sea continuum, the unstructured grid also includes the lower part of the major rivers flowing into the Adriatic Sea. The numerical grid consists of approximately 110 000 triangular elements with a resolution that varies from 5 km in the open sea to a few hundred metres along the coast and tens of metres in the inner lagoon channels (Ferrarin et al., 2019).

The SHYFEM simulations are forced by atmospheric fields from the global ECMWF reanalysis ERA5 (Hersbach et al., 2020), sea boundary and initial conditions from the Mediterranean Forecast System reanalysis (Tonani et al., 2009) and observed or climatological water discharges at the main river mouths.

4 Application examples

The ChemicalDrift functionalities are demonstrated through a few modelling examples described in the following. The presented simulations are not meant to provide conclusive quantitative results as some of the input parameters are uncertain and the model remains to be validated. The examples are run with OpenDrift 1.9.0 (e.g.) and may not be functional in the far future.

4.1 Organic pollutant discharge in Kattegat

A simple example of a ChemicalDrift simulation is presented with a description of the running script. The example is also available in the gallery section on the OpenDrift reference page (opendrift.github.io), where it is continuously updated with live forcing data. The simulation models the fate of a mass of an organic pollutant (phenanthrene) released outside the northern coast of Denmark over a period of 2 d. Simulation set-up is done as shown in the listing below by loading the OpenDrift modules, defining a ChemicalDrift instance, adding the necessary readers for forcing data and configuring a set of parameters. In detail, metocean forcing data including surface winds, ocean currents, sea water temperature and salinity are provided by the Norkyst800 model, while mixed-layer depth is set to a constant value of 40 m. Vertical mixing is activated, and the model used for the diffusivity profile is selected. Released chemicals are assumed to be 90 % in the dissolved fraction, and 10 % are adsorbed to particles. Degradation rates for phenanthrene are overridden in order to produce a clear effect in this short demo, and half-life constants are set to 6 h in the water column and 12 h in the sediment layer.

The model is configured at this stage. The next step is to seed the Lagrangian chemical elements and run the simulation. In this example, 500 chemical elements with a mass of 2000 µg are seeded within a radius of 2000 m around the selected position and within the upper 10 m of the water column.

Simulation results are illustrated in Fig. 2, showing how the chemical is partitioned between the dissolved, adsorbed-to-SPM and sediment compartments and advected in the northward direction. As expected, dissolved chemicals are mainly close to the surface and transported horizontally at higher velocities, while it can be seen that elements adsorbed to SPM travel for shorter distances at a higher depths. The distribution of mass between the three modelled compartments is shown in Fig. 3, where the exponential decay due to the combined effect of degradation and volatilization can also be observed. In detail, due to the short half-lives used in the simulation, more than 40 % of the initial mass was removed during the simulated period. Two clear strong resuspension events are illustrated after approximately 8 and 34 h, as shown by the rapid increase in the mass of chemicals associated with SPM at the expense of the sediment compartment followed by periods with a pre-dominant exchange in the opposite direction as elements representing suspended particles are sedimented again.

Figure 2 Simulated transport and fate of a mass of phenanthrene released outside the northern coast of Denmark, showing horizontal advection (a, b, c) and vertical distribution in the water column and at the sea floor (d, e, f). The colour of the Lagrangian elements indicates the partitioning between dissolved, adsorbed to SPM, and sediments.

Figure 3 Exponential decay and mass distribution across the dissolved, adsorbed-to-SPM, and sediment phases for phenanthrene.

4.2 Naphthalene and benzo(a)pyrene discharges from open-loop scrubbers

ChemicalDrift is demonstrated using input data from the STEAM model to simulate emissions of selected PAHs (i.e. naphthalene and benzo(a)pyrene) from open-loop scrubbers. The simulated region includes emissions for January 2019 in the area between 8 and 15^∘ E and between 53 and 60^∘ N. The simulation also covers an additional period of 2 months after the emission of the target chemicals ended, showing how the different chemical properties (Table 2) strongly influence their environmental fate. Results are illustrated in Figs. 4 and 5. In both figures we see that naphthalene is predominantly in the dissolved phase, while benzo(a)pyrene, which has a 3 orders of magnitude larger K_OW, is mostly adsorbed to particles. Initial discharges along ship lanes are observed in Fig. 4a and d. In the figure the diameter of each Lagrangian element is reduced as its associated mass decreases, indicating the target chemical's removal by degradation and volatilization. This is clearly observed for naphthalene (top right), while benzo(a)pyrene (bottom right) is almost entirely preserved in the simulated period since it has much longer half-life values and is to a large degree adsorbed to particles and more rapidly deposited to the sediment layer, where its degradation is much slower. The same conclusions about the different behaviours of the two chemicals are also observed in Fig. 5, which illustrates the two chemicals' vertical profiles in the water column and sediments (top) and the time evolution and phase partitioning of the mass of the chemicals (bottom).

Figure 4 Emissions of naphthalene (a, b, c) and benzo(a)pyrene (d, e, f) from open-loop scrubbers derived from the STEAM model for January 2019 in the sea region around Denmark. Element diameter is directly proportional to the associated mass. Transport and fate are calculated by ChemicalDrift over a 3-month period. Naphthalene is predominantly in the dissolved phase and more easily degraded and volatilized. Benzo(a)pyrene is mostly adsorbed to suspended particles and deposited to the sediment layer.

Figure 5 Emissions of naphthalene (a, c) and benzo(a)pyrene (b, d) from open-loop scrubbers derived from the STEAM model for January 2019 in the sea region around Denmark. Northward view of the distribution of chemicals in the water column and sea floor after 3 months (a, b) and evolution of the mass of chemicals in each model phase (c, d).

4.3 Multiple emission sources in the Adriatic Sea

Modelling of emissions from multiple sources is demonstrated. The target chemical considered is benzo(a)pyrene, and emissions include discharges estimated in August 2019 from open-loop scrubber data derived from STEAM and discharges from 49 rivers obtained from water-quality monitoring data and river loads.

Snapshots of the simulation are shown in Fig. 6, where colours are used to identify model compartments (panels a and b) and emission sources (i.e. ships or rivers, panels c and d). Emissions from scrubbers are to a large degree confined to the proximity of shipping lanes even at the end of the simulated period, and this is explained by also considering that benzo(a)pyrene has a high K_OW (Table 2) and is strongly associated with the solid fraction which is rapidly deposited to the sediment layer in proximity to the ship tracks. For the same reason, emissions from the rivers are mostly located along the shallow coastal regions and slowly transported southwards by the cyclonic Adriatic gyres.

Simulation results are further illustrated in Fig. 7, where the mass of target chemicals and the corresponding mass removed by degradation and volatilization are plotted versus time. The total emitted mass of chemicals in August 2019 is 17 500 g, and only about 500 g is removed by the end of September 2019. Degradation is dominant (470 g) compared to volatilization (30 g), which is also explained by the strong affinity of benzo(a)pyrene with the solid fraction, which does not volatilize. Most of the degraded mass is removed in the dissolved phase, even if most of the chemicals are associated with the sediment layer, which is explained since $t_{\mathrm{1}/\mathrm{2}\mathrm{W}\mathrm{tot}}$ is approximately 30 times shorter than $t_{\mathrm{1}/\mathrm{2}\mathrm{S}\mathrm{tot}}$ (Table 2).

The example is also utilized to illustrate the dependency on temperature and salinity of the solid–water partition coefficient which is calculated according to Eqs. (8)–(10). This is reflected e.g. in the desorption rate from the sediment to the dissolved phase, k₄₁ (Eq. 2), as shown in Fig. 8.

Figure 6 Emissions of benzo(a)pyrene from both open-loop scrubbers and rivers in the northern Adriatic Sea for August 2019, simulated in ChemicalDrift over a 2-month period, showing how the use of Lagrangian modelling allows for seamless separation of the different sources. Two time steps are shown, 8 August 2019 (a, c) and 29 September 2019 (b, d). The colours indicate the phase partitioning (a, b) and the pollutant source (c, d).

Figure 7 Emissions of benzo(a)pyrene from open-loop scrubbers and rivers in the northern Adriatic Sea. Total emitted mass (a) and removed mass by degradation and volatilization (b). Mass in water column refers to all phases in the water column, including the dissolved and SPM phases.

Figure 8 Effect on the desorption rate k₄₁ of benzo(a)pyrene from the sediments to the dissolved phase, due to SalinityCorr (Eq. 9 a) and TempCorr (Eq. 10 b), produced by the simulation of multiple source emissions in the northern Adriatic Sea.

4.4 Modelling of open-loop scrubber emissions in European seas

ChemicalDrift is applied to open-loop scrubber emissions for the entire European region from January to March 2019. The target chemical is phenanthrene, and fate and transport modelling are calculated until December 2019. This example is utilized to demonstrate the method write_netcdf_chemical_density_map(…), which is implemented in ChemicalDrift by adapting the corresponding method in the Radionuclides module. The method calculates gridded concentrations of the target chemical and exports these to a netCDF file. Gridding is applied as a post-processing step at the end of the simulation, providing average concentrations (10×10 km, user-definable) of the target chemical over the whole domain for each of the modelled phases. Estimated environmental concentrations are shown in Fig. 9. The time step for the averaging in each grid box is also user-definable. For the modelled phases in the water column, e.g. the dissolved and adsorbed-to-SPM phases, it is possible to define vertical grid size or to consider as in this case the whole column from the surface to the sea floor. Concentrations are saved to netCDF files. Results for the sediment phase are plotted (µg kg⁻¹), where the mass of the chemical is the mass of Lagrangian elements deposited on the sediment layer in each grid box, and the total solid mass is calculated from user-definable sediment density and active-layer depth and porosity.

We see that only a few ships with scrubbers were operating in the Mediterranean in 2019, since it was before global sulfur cap regulation. Dissolved chemicals are more diffused, while chemicals attached to particles sink and have smaller lateral diffusion, and hence ship routes are more easily observed.

Figure 9 Mean concentrations of phenanthrene emitted from open-loop scrubbers (January–March 2019). Concentrations are calculated separately for the dissolved (a), adsorbed to SPM (b), and sediments fractions (c).

5 Conclusions

ChemicalDrift, a new Lagrangian model for transport and fate of chemicals in the aquatic environment, has been presented. The model is implemented as a new module and is fully integrated within the open-source framework OpenDrift in order to combine the newly implemented chemical processes with the framework's advanced hydrodynamical capabilities and to provide a flexible interface with most of the available metocean input data sources.

The modelled chemical processes include a dynamic partitioning between the different phases that pollutants can be associated with in the aquatic environment, chemical removal by degradation and volatilization, and sedimentation and resuspension of chemicals associated with suspended particles. Target chemicals are modelled as Lagrangian elements that are transported and exposed to changing environmental conditions, e.g. temperature and salinity. The dependencies of chemical processes on temperature and salinity changes are formulated and implemented.

The focus of the presented work has been on modelling organic pollutants in the marine and coastal regions. The model functionalities are demonstrated through a sequence of simulation examples. The presented examples are only for demonstration purposes, providing insights into the combined effect of the modelled physical and chemical processes and presenting the potential of the proposed model. Accurate tuning of input parameters, sensitivity analysis, model validation and quantitative modelling results are deferred to future publications. Specifically, in the scope of the Horizion 2020 EMERGE project, ChemicalDrift is planned for use in calculating concentrations of different pollutants, including PAHs and heavy metals, both at the European regional scale and at a finer scale in the northern Adriatic Sea and in the Øresund strait. Additionally, the model can also be applied to other case-study areas in the EMERGE project, namely the Piraeus port-case study and the Aveiro Lagoon case study, where the Delft3D suite will be utilized independently by other project partners, allowing one to compare the results between the two models. Datasets with measured concentrations of organic pollutants in sediments obtained from monitoring programmes in the Baltic region and in the northern Adriatic are identified and available and can be used in the future to attempt validation of the model.

Additional functionalities are also planned for the future, including support for hydrolysis, photolysis and biodegradation as distinct sub-processes, support for non-ionizable organic compounds and a more advanced resuspension scheme which will include the contribution of wave-induced stress to resuspension.

High flexibility is demonstrated by the presented examples, utilizing a selection of different input sources and at different spatio-temporal scales: an interface to the STEAM model for shipping emissions, including support for open-loop and close-loop scrubbers, is implemented and demonstrated. The Lagrangian framework offers seamless tracking and separation of chemicals emitted from different sources such as emissions from shipping or from rivers. The modelled chemical processes depend on a relatively large set of parameters, including e.g. the solid–water partitioning coefficient, the Henry law coefficients and the overall degradation rates. A database of chemical parameters for a set of PAHs is integrated into ChemicalDrift, and values for the three PAHs are presented.