Combined integral and particle model for describing the dispersion, dilution, terminal layer formation and influence area from a point source discharge into a water body

Contamination of coastal water is a persistent threat to ecosystems around the world. In this study, a novel model for describing the dispersion, dilution, terminal layer formation and influence area from a point source discharge into a water body is presented and compared with field measured data. The model is a Combined Integral and Particle model (CIPMO). In the initial stage, the motion, dispersion and dilution of a buoyant jet are calculated. The output from the buoyant jet model is then coupled with a Lagrangian Advection and Diffusion model describing the far-field. CIPMO ensures that both the near- and far-field processes are adequately resolved. The model either uses empirical data or collects environmental forcing data from open source hydrodynamic models with high spatial and temporal resolution. The method for coupling the near-field buoyant jet and the particle tracking model is described and the output is discussed. The model shows good results when compared with measurements from a field study. A combined integral and particle model (CIPMO) has been designed to solve transition zones of plume jet trajectories by coupling the near-field and far-field with a transition function The developed model system can use measured data and hydrodynamic modelled data to predict the terminal layer formation, spread of the particles and dilution from point source discharges A comparative analysis of a tracer tracking experiment demonstrates that CIPMO simulations coincide well with the in situ observations of transport and dilution CIPMO can be applied as a decision support system for environmental impact assessments of discharges from industrial operations, landfills, and wastewater treatment plants A combined integral and particle model (CIPMO) has been designed to solve transition zones of plume jet trajectories by coupling the near-field and far-field with a transition function The developed model system can use measured data and hydrodynamic modelled data to predict the terminal layer formation, spread of the particles and dilution from point source discharges A comparative analysis of a tracer tracking experiment demonstrates that CIPMO simulations coincide well with the in situ observations of transport and dilution CIPMO can be applied as a decision support system for environmental impact assessments of discharges from industrial operations, landfills, and wastewater treatment plants

1 Introduction into water, and how the results can be used to understand how a discharge will affect the receiving waters, including its flora and fauna. The model output is compared to an in situ study conducted on Langøya, in the outer Oslofjord in 2008 (hereafter referred to as the Langøya study) [14], and the results are discussed.

Development of CIPMO
CIPMO incorporates four models that are designed to solve transition zones of plume jet trajectories. A schema describing the different steps are presented in Fig. 1.
The first model, denoted the near-field model (A) consist of an initialization model i.e. the Zone of Flow Establishment (ZOFE) [15], denoted (A. 1), and an integral model [16] (A.2). The ZOFE model (A.1) describes the region adjacent to the discharge [17], which determine the initial conditions of the sequential model. The integral model (A.2) computes the trajectory and concentration of a buoyant jet discharged into a stable density stratified environment, i.e. coastal waters and lakes.
The near-field model (A) describes the near-zone of the buoyant jet (see Fig.1), and is coupled to the Lagrangian particle tracking model LADIM (C) 2 [18], with the coupling model denoted (B). The lagrangian particle tracking model describes the far-field zone (see Fig. 1).
The particle tracking model LADIM (C) [18], developed by the Institute of Marine Research in Norway, uses a numerical technique to estimate the convective and diffusive motion from a large number of fluid elements that represents a continuous fluid or dispersed particles.
The last model, the Kernel Density Estimation Model (D), estimates the probability density function of the dispersed particles to quantify the impact of the discharge.
CIPMO is programmed in Python which is an interpreted, high-level, general-purpose programming language. In order to solve the challenge of a dynamic environment, CIPMO collects environmental data, like salinity, temperature and currents, from hydrodynamic models with high spatial and temporal resolution. Forcing data can be either empirical or retrieved from the Regional Ocean Modeling System (ROMS). ROMS is a free-surface, terrain-following, primitive equations ocean model, [19,20].
The output of CIPMO can be visualized as animations, images, or as graphs. The graphics can represent either transient snapshots of spesific time-steps, or as concatenated data from multiple simulations. The concatenated illustrations are accessible as either probabilistic or as arithmetic mean, however the raw-data are accessible for multivariate statistical analysis. Hence, a multitude of statistical representations in graph-form e.g. maximum, minimum, percentile and correlations are readily available.

The near-field model (A)
The integral near-field model (A.2) in CIPMO is similar to the CorJet model, 3 differing with its implementation of the ZOFE model (A.1) [15]. Due to the assumption of axisymmetric self-similarity it is necessary to implement a ZOFE, as it dictates the initial The near-field model (a) describes the near-zone of the buoyant jet and is coupled to the Lagrangian particle tracking model LADIM (c) which describes the far-field, with the coupling model denoted (b). The last model, the Kernel Density Estimation Model (d) quantifies the impact of the discharge trajectory of the buoyant jet plume until the buoyant plume satisfies the axisymmetric selfsimilarity condition, and the integral model is valid.
The integral model (A.2) consists of 10 coupled ordinary differential equations formulated to describe mass, momentum, buoyancy and scalar quantities in the turbulent jet flow [16]. This system of ordinary differential equations are implemented in the CorJet. In CIPMO, the system of first order differential equations are solved numerically with the Scipy 4 implementation of the 4th order Runge-Kutta solver. The integral model uses an entrainment closure that differentiates the different contributions of azimuthal shear mechanisms, e.g. advected momentum puff, and transverse shear, e.g. wakes, plumes and jets. The numerical solution of the system of ordinary differential equations predicts the averaged dilution, geometrical position of the plume trajectory and dispersion of plume width at the terminal layer, i.e. the instance the density gradient of the buoyant plume and the stratified environment becomes infinitesimally small. At this instance the solution exhibits strong instability, i.e. stiffness issues and the integral model is no longer valid to determine the proceeding dilution and plume trajectory.

Coupling between the near-field and far-field models (B)
To assert continuity between the near-field and the far-field models the initial far-field particle position at the transition point are distributed with a transition function, i.e. a coupling. The main scope of the transition function is to ensure that the far-field model conserves the center position and the degree of dilution at the transition.
The coupling model consists of four independent models: (B.1) a modified jet model based on Jirka [21], (B.2) a horizontal diffusivity model based on Okubo [22], (B.3), a vertical diffusivity model based on Forryan [23], and (B.4) a Monte Carlo particle allocation model. The three first models address and resolve excess momentum from the near-field model, conserves the continuity of mass, and incorporates both horizontal and vertical diffusivity. They are then integrated together using the Monte Carlo model. Since the coupling model singularly transfers data from the near-field to the far-field it must be categorized as a one-way coupling model.

Jet model (B.1)
At the termination point of the near-field model, the buoyant jet can still possess substantial momentum, even though the termination criteria are satisfied. It is imperative that the surplus momentum from the jet is conserved to avoid an unrealistically conservative evaluation of dilution. Hence, a modified 2-dimensional integral jet model based on Jirka [21] is implemented into the coupling model to conserve energy, momentum and continuity of the jet plume from the near-field model.
The shape of the jet is in accordance with the Laputz et al. [24] methodology retransformed from the circular shape determined by the effluent pipe evaluated throughout the near-field and into a rectangular form. This reshaping is consistent with the magnitude ratio between vertical and horizontal diffusivity from the ambient water and the vertical entrainment resistance since the jet is in density equilibrium with the surrounding water.

3
The modified model consists of five first order differential equations and a set of support functions. The set of differential equations are resolved with the Runge-Kutta method that is incorporated into the Scipy package ODEINT 5 [25].
The major modifications implemented in the coupling model consists of removing the buoyancy terms, since these terms are solved within the near-field model. This affects the continuity equation: where q is the volume flux, s is the arc length of the buoyant jets centre point trajectory, and e is the entrainment rate. In a similar manner to the near-field model, the jet model employs an entrainment closure scheme that differentiates between the respective contributions of transverse shear and of internal instability mechanisms: where u c is the excess velocity, 1s is the entrainment coefficient for a pure jet, 3s the pure wake, 4s is the entrainment coefficient for a pure advected plane puff, u a is the ambient velocity and is the horizontal angle which relates the x-axis to the centre trajectory of the buoyant plume.
In Eq. 2 the entrainment contributions from the pure plume, the thermal and the vertical angle ( 2s , 5s and , respectively) are unimplemented due to the absence of buoyancy within the coupling model.
The horizontal momentum components are still consistent with Jirka [21], however, the vertical angle is not present: where e is the entrainment function, u a is the ambient advective velocity, a is the angle between the ambient velocity vector and the x-axis. F d is the quadratic drag force, and i x and the i y are unit vectors which are further described in [16].
The positional equations are simplified, where the vertical angle is neglected, in addition to the Cartesian z-component: where x and y are the Cartesian positional coordinates, and account for the East-West and North-South directions respectively.
(2) e = 2u c ( 1s + 3s u a cos u c + u a + 4s u a sin |cos ) ds =eu a cos a + |F d |i y The width of the jet, b, is evaluated based on the integral variable: where M is the resultant momentum, q is the integrated continuity of the volume flux and b is the width of the buoyant jet.

Horizontal and vertical diffusivity (B.2 and B.3)
To avoid underestimating the entrainment contribution from the ambient water due to passive diffusion -two diffusion models, i.e. horizontal and vertical, are implemented within the coupling model. To account for horizontal diffusivity, the method described in Okubo [22] and in Matsuzaki and Fujita [26] is implemented. The methodology consist of curve fitting in situ data, from submerged dye and drifting buoy oceanic experiments, to power functions. The function implemented in the coupling model is consistent with the curve fitted function of Okubo [22]: where b is the jet width, of the dispersed medium, evaluated in Eq. 7 and K h is the horizontal diffusivity coefficient, K h .
The horizontal turbulent diffusion velocity, U h_diff , is calculated as follows: where t is the time-step, and R is a stochastic number in the range of [- To account for horizontal diffusivity the vertical diffusion model, described in Forryan et al. [23] is used. The model is formulated as a non-linear function, with a set of constants, which are fitted to in situ observations. The in situ data are sourced from three separate ocean regions in the North Atlantic and Southern Oceans. The model is optimized with regards to Richardson numbers greater than 1. It correlates significantly with observed mixing in the presence of mesoscale ocean features, as described by [23], and is formulated as the following: where s and s are constants, K bs is the background diffusivity, K os is the turbulent diffusion under neutral stability, and Ri is the Richardson number-which is a function of where N is the buoyancy frequency, and S h is the vertical shear.
The vertical turbulent diffusion velocity, U v d iff , is calculated as follows:

Monte-Carlo model (B.4)
The coupling model utilizes a Monte Carlo model to combine the respective sub-models, that is the jet model and the horizontal and vertical diffusivity models, into an merged model system. The Monte Carlo model allocates n-number of particles, which are free to propagate in 3-dimensions, at each time step. The n-number of particles are concatenated into evenly spaced numbers over a specified interval, i.e. a user input interface that determines the duration of the time-step, t .
The single determinator for the particle grouping mechanism is the o-number of jet length step o× s , that occurs within the specified time-step, t . The mean trajectory of the particles is determined by the spatial coordinates x, y as seen in Eqs. 5 and 6. z is the depth evaluated at the endpoint of the near-field model.
At each numerical length step, s , n-numbers of particles are displaced iteratively from the centre position according the following equations: R is an independent stochastic number which draws random numbers from a discrete uniform distribution, and s is the step length.
This generates a field of particles clustered within the spatial boundary of the jet and the diffusivity components from the ambient water. This coupling process is comparative to a process described by Nekouee et al. [1], where a near-field model and a far-field particle tracking model has been coupled. Although the near-field models are not similar, the techniques of Nekouee et al. [1] and CIPMO are comparable. This transition from a continuous jet to discrete particles enables the possibility to follow a representative selection of the buoyant jet's particles during their dispersal in the water body.
The number of particles is a multivariate function dependent on computer resources i.e., the number of CPUs, amount of RAM available, and the total quantity of time-steps.

LADIM (C)
LADIM [18] is the latest generation of the TRACE [1] subroutine originally coded by Jarle Berntsen (Institute of Marine Research, unpublished 1991) that approximates the trajectory and dispersion of particles from advection and diffusion schemes in three directions. The model represents the domain as an array of square grid cells. The horizontal advective velocities are sourced from the ROMS hydrodynamical model's current predictions, and interpolated from the grid centers to the grid square corners. The diffusivity velocities are parameterized with scaled stochastical components: where R is an independent stochastic number consistent with the description of Eq. 13, K h is the horizontal diffusivity coefficients and K v is the vertical diffusivity coefficients.
The velocity components are evaluated as following: The free surface is handled as a zero-flux boundary. Particles are constrained from permeating this boundary.

Boundary conditions
Land boundaries are treated, in addition to the free surface, as non-flux. However, the scheme developed by Ådlandsvik [18] 6 is implemented. This scheme is treating stationary particles by arbitrary reallocating the particles in close proximity to the stationary position and re-introduces the particles into the active environment. This reduces the erroneous numerical derived density build-up of particles nearshore.
Bottom boundaries are treated as non-flux boundaries, which reflects particles that are interacting with the bathymetric boundary by: where z s,i is the vertical position of a stationary particle and H s,i is the bathymetric depth evaluated at the same spatial position [18].
CIPMO has integrated an Individual Base Model (IBM) that is a collection of two alternating horizontal diffusion models with a dampening wall function, a vertical diffusivity model, and a scheme that accounts for tidal effects.

Horizontal diffusion scheme
The horizontal diffusions models are the Okubo model [22], which is described in chapter 2.2.2, and the horizontal Smagorinsky-Lilly turbulence model [27,28]. The Smagorinsky turbulence model approximates the turbulent diffusion based on the shear velocity components from adjacent cells: where C Smag is the Smagorinsky constant, x, y, u, v are, respectively, the decomposed horizontal positions and velocity components. The constant C Smag is adjusted to improve results and is usually evaluated between 0.1 and 0.2.
The Okubo model [22,26] is implemented in the IBM similarly as in the Coupling model (B): where L is calculated as the standard deviation of the positions of particles (defined as 3 ), based on Okubo [22]: where n p is a clustered quantity of particles analogous to the grouping order as seen in the coupling model. Hence, particles are designated a clustered group throughout the entire simulation period. This grouping scheme increases the dynamic behaviour and ensures that previous stochastic interaction for a specific cluster group is conserved in a continuous manner.
A brute force optimization scheme is implemented to evaluate the quantity of particles designated in each group. The main scope is to conserve the grouping identity evaluated in the coupling model. Since the length of each coupling model trajectory is statistically likely to diverge, the grouping order is not homogeneous for each time-step. Hence, a mean trajectory will be used as a proxy to concatenate the quantity of the particles for each group.
The quantity of particles in each group are evaluated based on the following optimization scheme: where n is the net accumulation of particles released at time-step t i+1 =t i + t , v is the mean advective ambient velocity evaluated at the terminal layer depth, t is the numerical timestep, s is the length step applied to solve the system of differential equations from the jet model, and is an integer value optimized to determine the group quantity.
Due to the high diversity of cell sizes applied in compatible driving force data, the Smagorinsky diffusion model's applicability for the initial evaluation of the horizontal diffusion in the far-field is suboptimal. It is imperative that the scale of the turbulent eddies is comparative to the spatial boundary of the passive effluent, otherwise the impact of the eddies will be recategorized into advective terms. The Smagorinsky coefficient C Smag , is applied to address potential discrepancy of cell sizes from the forcing data and scales of turbulent eddies with diffusive properties in relation to the effluent spatial extension. However, if a semi optimized value of the coefficient C Smag , is implemented throughout the simulation period, the accuracy of the scaling factor will diverge as the spatial extension of the effluent contracts or expands from a critical range. Minimize: Subject to: The following alternating function is implemented to mitigate some of these issues: This ensures that the most suitable horizontal diffusivity model is selected by CIPMO.

Dampening function
Since the Smagorinsky model determines the horizontal eddy diffusivity based on shear velocity intensity, the maximum value will usually occur at near solid boundaries. With this theoretical approach, unlike observational studies [29], turbulent eddies are not damped near the wall boundaries. Hence, a damping mechanism [28] is implemented to mitigate this erroneous behaviour: with y + evaluated as a function of where r is the distance to the wall, u is the friction velocities, and is the kinematic viscosity.

Vertical diffusion scheme
The vertical diffusion scheme is described in chapter 2.2.2 and is implemented into the IBM model identically.

Kernel density estimation model (D)
Kernel density estimation (KDE) [30] is used to calculate the concentration, vis-à-vis the dilution in the far-field of the allocated particles. KDE is a hybrid combination of discrete binning, i.e. grouping adjacent particles together within a discrete domain with an automatic locally adaptive smoothing kernel. The general principle of the concentration estimation method [31] used in CIPMO, is that n-particles of unifor mass are released, and each particle have uniform mass. The particles positions are given by {X 1 ,...,X n }, {Y 1 ,...,Y n } and {Z 1 ,...,Z n }. The concentration is evaluated by constructing a grid, with d-dimensions with the node position The resolution of the grid is dependent on the following optimization scheme: where x, y, z is the respective cells dimension in the three Cartesian directions, G is a Gaussian filter function, C threshold is either the predicted no-effect concentration (PNEC), EQS value or an analogous threshold value, is a safety factor evaluated as a constant between 1000 and 10,000, C p,i is the concentration evaluated for a single particle, and m p and n p are respectively the total accumulation of mass and particles released at an arbitrary time-step. The scope of the scheme is to allocate a singular particle, i, within a grid element with the undetermined dimensions x , y , z , where the particle represents a finite quantity of the total dispersed mass, i.e: By altering the proportional relationship between the total numbers of particles, n p , and the spatial dimensions x , y , z , the concentration evaluated within the grid-cell can be altered.
The optimization scheme utilizes this relationship to find the smallest grid cell which corresponds to a concentration of a singular particle which must be equal or smaller than the threshold value of the target chemical substance.
This mass grid dimensional discretization scheme ensures that the concentration field is adequately resolved. The concentration field will be artificially high if this scheme is not resolved properly.
Equation 26 has included the G term, a Gaussian filter function, which contributes to the grid element, x , y , z , as a smoothing agent.
The concentration of the particles is given by: where ρ(x, y, z) is the particle density function, and x y z is the volumetric content of the fluid per unit volume. The particle density function is evaluated on a multi-dimensional grid. The grid is based on the free software library NUMPY [32]. The particle density function ρ(x) is approximated via the KDE function: where K is the kernel, a positive function controlled by the bandwidth parameter h and X i , Y i , and Z i is the particle position. CIPMO uses the free software machine learning library (26) Minimize: Subject to: Scikit-learn sub-module "sklearn.neighbors.KernelDensity" [25] to evaluate the particle density, ρ(x), of n particles. CIPMO quantifies the potential environmental impact from the effluent sources in two distinct ways with an arithmetic and an probabilistic scheme. Both are based on evaluating the probability density function on a three-dimensional grid. They are in turn projected into a horizontal plane by selecting the highest numerical values from the vertical z-axes. However, the succeeding steps diverge.
The arithmetic scheme consists of a for-loop mechanism which iteratively summarizes all the projected horizontal meshes from a discretized selection of time-steps and/or independent simulations, and lastly divides the entire multiple with the quantity of iterated steps, i.e. evaluating the arithmetic mean of projected planes.
The probabilistic scheme is similarly formulated on a for-loop mechanism, but uses Boolean algebra. The numerical values from the horizontal plane are either denoted with the value zero or one, based on the following logical test: if the values are greater or equal to a threshold value, such as EQS or PNEC, it is assigned a value of 1, otherwise 0 is assigned. The looping mechanism iteratively summarizes all the planes designated from a set of time-steps and/or independent simulations. Lastly, the multiple is divided by the quantity of iterated steps. The output provides a result which can be interpreted as the probabilistic chance that a specific area will exceed a certain threshold value.

Model evaluation-comparison with field data
To investigate the predictive capabilities of CIPMO, the model is evaluated against a field experiment where discharge dispersion from Langøya landfill was studied [14]. Langøya landfill is operated by NOAH (Norsk avfallshåndtering), located on an island in the outer part of the Oslofjord (Fig. 2).
The Langøya study was performed by the Norwegian Institute of Water Research (NIVA) in February 2008. A tracer (Fluorescein, MS-200) was added over a period of 40 h, 100 l all together. The tracer tracking was done February 18 and 20, starting approximately 24 h after the tracer was first added. This resulted in two different measurement periods, period 1 lasting from 24 h-32 h, and period 2 lasting from 48 h-54 h, after the tracer was first added. 40 samples were taken during period 1 and 30 samples during period 2 (see Fig. 5).
The approach for the comparison study employs probabilistic modelling (i.e. simulating the release multiple times). This makes it possible to evaluate the influence of variable hydrological conditions such as currents, salinities, and temperatures during the release.

Input data and model setup
The pipe diameter, position and depth, the outlet flow velocity, and initial tracer concentration were replicated in CIPMO to be identical to the Langøya study (see Table 2). The hydrodynamic data (currents, salinity and temperature) measured during the experiment are point estimates and thus not suitable (and not publicly available) as input to modelling the release over a larger geographic area. CIPMO was configured with forcing data from the hydrodynamic model FjordOs, 7 [33,34]. Data from 14 randomly chosen time periods, distributed throughout February 2018 and 2020 (and one in 2019), were selected (Table 1). Each release is modelled for 54 h.
FjordOs is co-funded by nine cross-sectoral institutions and covers the Oslofjord in the South-East of Norway, with a variable horizontal resolution between 50 and 300 m and 42 terrain following vertical S-layers. The temporal resolution is one hour.
An illustration of the FjordOs dataset is presented in Fig. 3. The map shows a snapshot of the current at 30 m depth February 22 2020 (Fig. 3a), the current roses shows the Fig. 2 Illustration of the analysis area for the Langøya study. The release point is marked with a red circle and the location used to compare temperature and salinity data is marked with a blue circle average percentage of current direction toward each of the principal compass points for February 2018 and 2020 (n = 69 days) (Fig. 3b) and the selected dates (n = 14 ) (Fig. 3c). The prevailing current directions (north-westerly) and velocity (mainly below 6 cm/s) is similar, suggesting that the selected data is representative, at least for this area.
Water density is important input data to CIPMO. A comparison of salinity and temperature data measured February 15 at a station called Breiangen Vest (see Fig. 2) 8 with data from FjordOs is given in Fig. 4. Although the temperature, and partly salinity, is greatly influenced by weather conditions and local circular patterns, the modelled data from  Table 1). The data are taken from a measuring point South-East of Langøya at 30 m depth, marked with a blue circle FjordOs are comparable to the measured data. Both depict a relatively strong pycnocline, i.e. where the density gradient is greatest within a body of water, at around 40-50 m depth.
The near-field model (A) configuration is in alignment with the study setup, as seen in Table 2, except for the time discrepancies between the in situ study and the available  Table 2 Overview of the input data and the discharge parameters used in the comparison study (see text for details). The input data is identical, while the hydrodynamic data are not. The model ID refers to four main models in CIPMO (see Fig. 1). The outlet of the pipe was located 80 meters from shore, at 38 meters depth The determination of the time-step is based on a cost-benefit assessment with respect to interactive relationship between the net quantities of particles and the length of the time-step. Based on a small amount of iterations between these two parameters, the time-step of the couple model and LADIM model was set to 15 min. LADIM (C) is additionally configured with a Smagorinsky coefficient, C Smag , of 0.1, corresponding to the recommendations in Niceno and Hanjalic [27], due to the mitigative diffusive environment that is associated with fjord systems.
The grid of the KDE model (D), denoted in the table as [ x, y, z] KDE , has a spatial domain of 13 km × 13 km × 60 m. The horizontal cell sizes is 10 m in East-West direction and 10 m in North-South direction, and the vertical cell size is 1 m. The centroid of the nested grid is allocated to the release location at a depth of 38 m. The grid consist of 101 400 000 cells. The cell sizes and domain is selected on the basis of Eq. 26. The kernel applied with the KDE method is consistent with the Gaussian function in the Python package and class "sklearn.neighbors.KernelDensity" [25].

Results from the Langøya comparison
The distribution and concentration measured in the Langøya study is depicted in Fig. 5. The figure shows the horizontal distribution of the discharge for the two periods, 24-32 h, left and 48-54 h right. The concentration is summed up for each depth measured and projected to the area, so the unit becomes g/m 2 . The lowest and highest detectable concentration in the Langøya study was 0.01 g/m 2 and 1.2 g/m 2 , i.e. the concentration in the coloured areas is between 0.01 g/m 2 and 1.2 g/m 2 . We have, based on Staalstrøm et al. [14], estimated the area with concentration > 0.01 g/m 2 to be 1.2 km 2 for the first period (24-32 h) and 2.2 km 2 for the second period (48-54 h). In the Langøya study they assumed that the concentrations were slightly overestimated and the area was slightly underestimated for the first time period (24-32 h) and visa versa for the second time period (48-54 h).
The transport and dispersion of the tracer release modelled with CIPMO is illustrated in Fig. 6. The maps shows the mean area with concentration above 0.01 g/m 2 calculated from all time-steps within the two periods, i.e. for period 1 the mean area of 448 maps (areas) (8h × 4 time-steps× 14 simulations) and for period 2 from 336 maps (6h × 4 time-steps × Fig. 6 Mean transport and dispersion of the tracer release simulated with CIPMO for period 1 (a) and period 2 (b). The concentration is depth-integrated means (see text) and is estimated from all time-steps within the two periods, i.e. n = 416 for period 1 and n = 312 for period 2 Fig. 7 Histogram of estimated area (km 2 ) with mean concentration above 0.01 g/m 2 for period 1 (a) (n = 448) and period 2 (b) (n = 336). The labels on the x-axis show the upper limit of the bin. The bin interval for period 1 is 0.25 km 2 and the bin interval for period 2 is 0.40 km 2 14 simulations). The mean area for period 1 is 1.00 km 2 [95% CI: 0.33, 2.00], and the mean area for period 2 is 2.01 km [95% CI: 0.81, 3.02].
A visual comparison shows a good match between CIPMO simulations and the Langøya study, although less particles are calculated on the southern side of the island with CIPMO, compared with the in situ measures in the last period. The total mean maximum concentration estimated with CIPMO was 2.38 g/m 2 (range 1.44-8.31 g/m 2 ) for period 1 and 0.24 g/ m 2 (0.04-0.46 g/m 2 ) for period 2.
A histogram of the estimated areas is shown in Fig. 7. For period 1, 67% of the areas estimated with CIPMO is smaller than the area estimated in the Langøya study (and 33% is larger) and for period 2, 55% of the areas are smaller (and 45% is larger). The differences in the areas estimated by CIPMO is primarily due to variation in hydrological conditions. Influence areas, showing the probability that the concentration in the water column exceeds 0.01 g/m 2 are illustrated in Fig. 8. The areas are calculated from all time-steps within the two periods, and thus the area will have a larger extension than a single release of the tracer. The areas overlap with the field data and the likelihood increase closer to the release point. The probability that CIPMO estimates that the tracer is transported outside the geographical boundaries measured in the field study is around 20% . This is expected since the influence areas are based on different forcing data and includes areas with low probabilities of exceeding the threshold of 0.01 g/m 2 . Using a higher threshold value will result in smaller influence areas (one would typically use an EQS or PNEC in a risk assessment study).
How and at what depth the substance of interest is trapped in the water column, also called the terminal layer formation, is of great significance. This can be used to see if the substance released will reach the surface or not, and whether or not the substance will mix with other discharges in the same area, or come into contact with areas of interest. In the Langøya study, the terminal layer formation was detected to lay between 33 and 49 meters depth [14].
The terminal layer depths calculated by the near-field model in CIPMO is presented in Fig. 9). The average depth (all years) is 35 meters [95 %CI: 30,40] (Fig. 9). The estimated depths have a normal distribution and are relativity stable (standard deviation is 2.2 meters) within and between the different years investigated. The results from the calculated terminal layer formation by CIPMO coincide well with the measured terminal layer formation in the Langøya study.

Discussion
The coupling model of CIPMO is a combination of multiple models, collectively conserving the continuity, momentum and diffusion terms and transforming continuous fluid into discretized particles. The coupling model includes the ambient diffusivity terms, usually associated with the intermediate and far-field zone [1][2][3]6].
The main benefit by including these terms is that the model is no longer constrained to operate with small timescales. The alternative approach is to significantly reduce the time-step, however, this will affect the overall simulation time. Hence, the reason for using an ensemble of models is to incorporate independent entrainment contributions which collectively increase the robustness, and consequently emulate an improved depiction of the actual phenomenon.
A drawback with the nearshore scheme implemented in CIPMO is the low success rate of repositioning stranded particles back into the ocean. Assumed a particle is stranded on a shore, the initial probability of successfully reposition of the particle is a function of cell dimension and particle position within the cell. Hence, if the particle is tangentially located at the land-water boundary it is likely that the displacement scheme is successful. However, if the particle is repositioned in a diametrical direction of the active environment, this mitigates a successful outcome in the next iterative step, and can potentially create a cascading loop of unsuccessful reallocation attempts. The positive aspect of the scheme is the numerical computation cost benefit this model has on the total simulation time of LADIM. An alternative approach is to implement a weighted function that utilizes a distance gradient from shore scheme. This can significantly increase the success rate of reallocating stranded particle back offshore, however, this may have a considerable numerical cost.
An alternative approach, as seen in [1], is to kill off the particles that are interacting with the boundary, instead of repositioning them. This approach will certainly reduce the quantity of stranded particles at stagnation points, such as the shoreline. This is a reasonable scheme if the stagnant condition of the particles is singularly a numerical artefact. However, it is difficult to justifiably retain this assumption if a water soluble component is considered, which is mainly the case with CIPMO.
When the fluid elements in LADIM interact with the bathymetric boundary, a reflective mechanism is activated and reallocates the element slightly above the bathymetric boundary. Based on the assumption that the fluid elements are in density equilibrium with the surrounding ambient water, it is important to assess whether this mechanism is too crude and should, as a minor improvement, be a function of the Richardson number.
An alternative pathway for fluid elements to interact with the jagged structure of the bathymetry is that the horizontal trajectory of these elements instead deflects, and that the vertical position remains unaltered. An alternating function which, based on the Richardson number, selects one of these two dichotomic mechanisms can alternatively be utilized. The Richardson number, Ri, is a dimensional number that expresses the ratio of the buoyancy term to the flow shear term. If the Ri < 0.25 the shear velocity can overcome the tendency of a stratified fluid to remain stratified, and hence the vertical deflective boundary mechanism can be deployed, otherwise the horizontal deflection is selected. It is a valid point to assume that upwelling can be interpreted as a causative attribute of the vertical deflection mechanism and should be included as a determining variable. However, the velocity component of upwelling is usually several scales of magnitude smaller than the horizontal velocity components and can minutely account for the system responds from Eq. 18.
The concentration field, depicted in Fig. 6, is evaluated by selecting the cells with the highest concentration value in Z-direction, as described in chapter 2.4. The end result is a two-dimensional array projecting the maximum concentration field in a horizontal plane. Hence, the result is intrinsically conservative, and depicts the most severe and perhaps a distorted concentration projection. It is important to realize that cells in proximity can diverge notably in altitude.
An alternative approach is to project either the median, mean or a percentile of a sorted group, from the vertical axes, however, this approach seems to counter a conservative approach, even though the graphical depiction can justifiably be termed somewhat more realistic.
The length of the numerical time-step, evaluated in LADIM and the coupling model, should at a bare minimum, not exceed the numerical time-step in the forcing data [35,36], which in FjordOs is 1 h. However, an adequate study to determine the causative effect that the duration of the numerical time-step has on either the concentration, or trajectory of the particles has not been conducted. Hence, the reasoning for applying the time-step, 15 minutes in the test case, cannot be properly justified with respect to a specific cause-and-effect analysis of LADIM and is based on results from a minor sensitivity study. The numerical time-step have multiple profound effects on the dynamics between the coupling model and LADIM. If the time-step is too small the excess conserved momentum in the jet sub model might be too pronounced. Consequently, a premature transitioning intersection to passive particles will result in underestimating the net entrainment. Since the coupling model does not include a comparative shoreline handling mechanism it will not re-allocate particles beached on the shoreline. To compensate for this it is necessary to constrain the magnitude of the time-step.
When compared with the Langøya study, there is little deviation between the measured values and CIPMOs calculations. The terminal layer formation, the spread of the particles and the dilution coincides well. A statistical modelling approach combined with hydrological models with high spatial and temporal resolution is advantageous when defining mixing zones (and other potential risk areas) and will give valuable input data to wastewater treatment planning as well as impact and risk assessments. Linking the output to an ecological model, or introducing so called valued ecosystem components in the modelling would further increase the usefulness of the CIPMO model. This study singularly utilizes modelled forcing data, however, CIPMO is compatible with multiple in situ and data assimilated models. In situ data of current and densitometric data in four dimensions, including a temporal dimension, are sparsely available, and never for large geographical areas. Hence, complex analysis of current trajectories is solely dependent on complex hydrodynamical models for the majority of coastal environments.
Hydrodynamic models provide continuous, detailed quantifications of ocean current, salinity and temperature. It is, accordingly, imperative that the modelled data are accurate. The FjordOs hydrodynamic data has been evaluated and found robust and accurate when compared to measured data in the Oslofjord [34].
The quantity of data enables CIPMO to provide detailed analysis, however, the accuracy of the hydrodynamical data is essential to quantify, otherwise the model will be prone to eventually cause erroneous conclusions. With good forcing data, such as FjordOs, CIPMO is capable to provide detailed analysis of releases into the environment.

Conclusion
This paper presents a description and evaluation of a model system that unifies the smallscale entrainment features from integral models with the predictive capabilities of Lagrangian particle tracking models. The coupling was accomplished by unifying four independent models, to account for conserved momentum at the termination point of the nearfield model, horizontal and vertical diffusion from the ambient water, and a Monte Carlo scheme to allocate the particles within the geometrical confinement of the continuous jet-diffusion system.
A central motivation is to utilize open-source output from hydrodynamical models as the main forcing data and solving the coupling and far-field computation offline to reduce the computational cost.
The comparative analysis of the empirical dye study at Langøya shows the ability of CIPMO to predict terminal layer formation, the spread of the particles and dilution from point source discharges. The study outcome demonstrates that the model system can be applied as a decision support system for environmental impact assessment from continuous, and abrupt discharge operations, derived from industries, landfills and wastewater treatment plants.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.