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 ﬂexible interface compatible with a wide range of available meto-cean data is made accessible by the integration, making the tool easily adaptable to different spatio-temporal scales and ﬁt 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 treat-ment plants. Speciﬁc 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 re-sults 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 ﬂexibility, 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.


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, andvolcanic 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., 2019Calgaro et al., , 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.

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 Open-Drift 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.

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.

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.2for estimating the adsorption and desorption kinetic constants, k ads (L kg −1 h −1 ) and k des (h −1 ).
Transfer rates from the dissolved phase to the DOC, SPM and sediment compartments (i.e.k 12 , k 13 and k 14 ) and the corresponding desorption rates (i.e.k 21 , k 31 and k 41 ) are given by Here the indices refer to the numbering illustrated in Fig. 1: C DOC is the concentration of DOC in the water column (kg L −1 ), C SPM is the concentration of SPM in the water column (kg L −1 ), ρ sed is the density of sediment solids (kg m −3 ), 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 sorptiondesorption 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 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). https://doi.org/10.5194/gmd-16-2477-2023 Geosci.Model Dev., 16, 2477-2494, 2023 Table 1.Parameters for the calculation of sediment-layer adsorption and desorption rates.

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 −1 ) is the ratio between the adsorption and desorption kinetic constants: (3) 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 2 = 0.87), following 1/k des ≈ 0.03K d .Using this approximation and the definition of K d , the adsorption and desorption kinetic constants can be estimated as (4) 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 −1 OC ) and the organic carbon content of the dry matter (f OC -g OC g −1 ).
Typical values for f OC are in the range 0.01-0.10g OC g −1 (Seiter et al., 2004).
Relationships between K OC and the octanol-water partition coefficient K OW (L kg −1 ) 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.K OC,DOC ≈ 2.88(K OW ) 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 −1 ), and C salt is the molar concentration of the salts in solution (mol L −1 ).Salt concentration is calculated starting from salinity considering the molecular weight of seawater salt, MolWt salt = 68.35g mol −1 (Schwarzenbach et al., 2016), and the water density, ρ.
The correction factor for temperature is based on the Arrhenius equation where H (J mol −1 ) is the enthalpy of adsorption for the given chemical and R is the universal gas constant (R = 8.3145 J mol −1 K −1 ).The model utilizes distinct enthalpies for DOC and SPM, indicated as H KOC,DOC and H KOC,SPM , respectively.

Sedimentation and resuspension
Particle sedimentation and resuspension dynamics in Chem-icalDrift 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 res + N (z 0,res,unc ), where z res is the resuspension depth, z res,unc is the resuspension depth uncertainty and N (0, z res,unc ) is a randomly generated number from a normal distribution with zero mean and standard deviation z res,unc (Simonsen et al., 2019a).

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 1/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 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.from the corresponding half-life constants t 1/2W tot , t 1/2S tot applying k = − log(0.5)/t1/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).

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 −3 s −1 ) can be calculated based on the properties of the two-film interface and the concentrations of the chemical: Here, C water (kg m −3 ) is the concentration of the chemical in the water phase, C sat (kg m −3 ) 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 −1 ) 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.∂ ∂t C water = k vol C water ) that is solved by C water (t 2 ) = C water (t 1 ) e −k vol (t 2 −t 1 ) , where the volatilization rate constant (s −1 ) can be calculated from Eq. ( 14).
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 sixring 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.
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 + t) = m i (t) e −k vol 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): Here, MTC a and MTC w (m s −1 ) are the air-side and waterside 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)  Here, MolWt (g mol −1 ) is the molecular weight of the target chemical, MolWt CO 2 is the molecular weight of CO 2 (44 g mol −1 ), MolWt H 2 O is the molecular weight of water (18 g mol −1 ), U 10 is the wind speed 10 m above the water surface (m s −1 ) and Sc H 2 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 Here, R is the universal gas constant (8.206 × 10 −5 atm m 3 mol −1 K −1 ), 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 3 mol −1 ), which is either specified by the user or estimated by where P sat (atm) is the vapour pressure of the target chemical at the reference temperature, C sat (g m −3 ) is the solubility of the target chemical at the reference temperature, and TempCorr(•, •, •) 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 .
3 Emissions, chemical parameters and metocean forcing

Ship Traffic Emission Assessment Model (STEAM)
Emission fields from STEAM by Jalkanen et al. (2021) 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).

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

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 hor-izontal velocities, mixed-layer depth, ocean surface winds, and SPM concentration, as summarized in Table 3. 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.

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.

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 mixedlayer 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 pro-duce 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, adsorbedto-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.charges 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 1/2W tot is approximately 30 times shorter than t 1/2S tot (Table 2).

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 postprocessing 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 −1 ), 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.

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 ca-   pabilities 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.

Figure 1 .
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.

Listing 2 .
Seeding elements and running simulation.

4. 2
Naphthalene and benzo(a)pyrene discharges from open-loop scrubbersChemicalDrift 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 (

Figure 2 .
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 .
Figure 3. Exponential decay and mass distribution across the dissolved, adsorbed-to-SPM, and sediment phases for phenanthrene.

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 dis-

Figure 4 .
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 .
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).

Figure 6 .
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 .Figure 8 .
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.

Table 1 .
Galí et al. (2022) 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 byGalí et al. (2022).

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

Table 3 .
CMEMS products and data variables.
Listing 1.Chemical simulation set-up and configuration.
Table 2) strongly influence their environmental fate.Results are illustrated in Figs.4and 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 https://doi.org/10.5194/gmd-16-2477-2023Geosci.Model Dev., 16, 2477-2494, 2023