Effects of Wave-Induced Processes in a Coupled Wave–Ocean Model on Particle Transport Simulations

Abstract

This study investigates the effects of wind–wave processes in a coupled wave–ocean circulation model on Lagrangian transport simulations. Drifters deployed in the southern North Sea from May to June 2015 are used. The Eulerian currents are obtained by simulation from the coupled circulation model (NEMO) and the wave model (WAM), as well as a stand-alone NEMO circulation model. The wave–current interaction processes are the momentum and energy sea state dependent fluxes, wave-induced mixing and Stokes–Coriolis forcing. The Lagrangian transport model sensitivity to these wave-induced processes in NEMO is quantified using a particle drift model. Wind waves act as a reservoir for energy and momentum. In the coupled wave–ocean circulation model, the momentum that is transferred into the ocean model is considered as a fraction of the total flux that goes directly to the currents plus the momentum lost from wave dissipation. Additional sensitivity studies are performed to assess the potential contribution of windage on the Lagrangian model performance. Wave-induced drift is found to significantly affect the particle transport in the upper ocean. The skill of particle transport simulations depends on wave–ocean circulation interaction processes. The model simulations were assessed using drifter and high-frequency (HF) radar observations. The analysis of the model reveals that Eulerian currents produced by introducing wave-induced parameterization into the ocean model are essential for improving particle transport simulations. The results show that coupled wave–circulation models may improve transport simulations of marine litter, oil spills, larval drift or transport of biological materials.

1. Introduction

A rapid increase in marine litter in the ocean has recently been recognized as a serious environmental problem. The role of the physical factors contributing to it (e.g., atmosphere–ocean–wave interaction processes) has been not yet fully understood. Lagrangian analyses represent the natural approach to studying oceanic transport based on model simulations and observational data [1,2]. The transport and accumulation of floating marine debris and the assessment of different scenarios of marine plastic distribution, along with a comprehensive synopsis of the tools currently available for tracking virtual particles. A detailed review of the Lagrangian ocean analysis, associated problems, sources of errors and validation issues was presented in [2,3]. Still, an improved understanding of the physical processes influencing the transport of particles is required [2]. The Lagrangian simulation can be assessed by performing a time-evolving analysis of the separation distance between the real track and the simulated ones. A skill score, based on the separation distance normalized by the length of the trajectory, has recently been proposed [4]. The separation distance between model simulations and observed trajectories has been estimated, showing that one day after initialization the distance was about 15–25 km, and five days later this increased to about 60–180 km. It was found that the root-mean-square error (RMSE) 13 days after the start of the integration was about 5 km [5]. A separation distance from model simulations versus observed trajectories in the first days after initialization of about 15 km was found [6].

Sea-state dependent processes affect the ocean circulation and thus also the results from Lagrangian transport models. A considerably enhanced momentum transfer from the atmosphere to the wave field is found [7] during growing sea state (young sea). Recently, a wind stress formulation depending on wind stress and wind–wave momentum released to the ocean was proposed [8]. Swell waves can even cause momentum transfer from the ocean to the atmosphere [9]. In growing sea states, waves extract momentum from the atmosphere, so that the ocean receives less momentum from the atmosphere than if waves were not considered [10]. Using stand-alone ocean or atmosphere models, the surface waves that represent the air–sea interface are not taken into account. This can cause biases in the upper ocean due to insufficient or, in some cases, too strong mixing [11], or because the momentum transfer is shifted in time and space compared to how the fluxes would behave in the presence of waves. Several parameterizations were recently proposed for momentum flux that is sea state-dependent, e.g., [12,13]. Recently, the role of the Stokes drift and wave-induced transport of floating marine litter was studied in [14], showing that accounting for the wave-induced Eulerian-mean flow significantly alters predictions of transport of floating marine litter by waves. The skill scores between the model and observations were improved by adding the Stokes drift in [15], postulating that for Lagrangian simulations, the Stokes drift’s contribution can be at the same order or even higher than the accuracy of the Eulerian circulation. However, the accuracy of Lagrangian simulations also depends on the hydrodynamic and Lagrangian model, demonstrating the need to tune Lagrangian models for the specific setup [16].

In a large and complex system such as the North Sea (Figure 1), minor perturbations may displace particles to very different drift regimes, causing strongly divergent particle trajectories [17]. Therefore, the particle distribution analysis is not intended to fully reproduce or explain the observations (the latter are extremely limited, particularly during extreme weather conditions). Instead, we aim to indicate the differences that can be due to wave-induced forcing processes. Although tide and wind-driven circulation in the North Sea seem to have been widely studied in the past, further research from a Lagrangian perspective with respect to coupling with waves and biogeochemical processes may be particularly relevant.

The importance of wave forcing for ocean circulation and sea-level predictions has been demonstrated [18,19,20,21,22,23], showing that the predictive skill of ocean circulation and sea level could be significantly enhanced by considering wave-induced processes. In extreme storm surge conditions over the North Sea, due to the strong non-linearity of wave–ocean–tidal interactions, wave–ocean coupling is considered to be significant for correct model predictions [18,24]. The impact of Stokes-related drift effects (Stokes–Coriolis forcing and Stokes drift advection on tracers and mass) were studied on the North Sea and Baltic Sea regions [25]. For tracer distribution and upwelling, the direct sea-state dependent momentum and energy fluxes are of higher importance than the Stokes drift processes implemented in the NEMO circulation model [23].

In this study, we investigate the role of wave-induced interaction processes in a fully coupled ocean–wind wave model system on Lagrangian transport modelling. Data from drifters deployed in the southern North Sea were used to assess the particle model simulations [26]. Lagrangian simulations were studied [27] based on the same drifter observations as in the present study and corresponding Lagrangian model simulations by using two circulation models. While only a stand-alone ocean model was used in the previous study [28], here, the same Ocean General Circulation Model (OGCM) is coupled to a wave model. The observational drifters [26] used to assess our simulations have previously been described in detail [27,29].

This paper is organized as follows. In Section 2, we describe the circulation, wave and Lagrangian transport models and the experimental setup. The evaluation of the model simulations is described in Section 3. Further, we assess the model results with direct comparisons of the simulations with the drifter data (Section 4). The drifter trajectories were modelled to investigate the importance of the wave effects, e.g., for search-and-rescue applications. We performed sensitivity experiments to investigate the impact of wave-induced forcing in the ocean model. The Lagrangian model used only the Eulerian current as provided by the hydrodynamical model simulations, by neglecting or considering the contributions from the Stokes drift or the wind drift corrections. The wave effects on the general circulation in the North Sea are studied. The paper ends with a discussion of the findings in Section 5 and the conclusions in Section 6.

2. Methods

2.1. Models and Set-Up

The Geesthacht coupled coastal model system (GCOAST) [30,31] was built upon a flexible and comprehensive coupled model system, integrating the most important key components of regional and coastal models. GCOAST encompasses (i) atmosphere-ocean–wave interactions, (ii) dynamics and fluxes in the land–sea transition, and (iii) coupling of the marine hydrosphere and biosphere. In our study, we used the GCOAST circulation, wave and drift model components to investigate the role of coupling in particle transport simulations in the North Sea. Those particles can be considered, for example, as simple representations of either oil fractions, fish larvae or search-and-rescue objects [32,33]. The wave–current interaction processes are momentum and energy sea state dependent fluxes, wave-induced mixing and Stokes–Coriolis forcing.

2.1.1. The Circulation Model NEMO

NEMO (Nucleus for European Modelling of the Ocean, [34]) is a framework of ocean-related computing engines, from which we use the OPA (Océan Parallélisé) package (for the ocean dynamics and thermodynamics) and the LIM3 (Louvain-la-Neuve Sea Ice Model) sea-ice dynamics and thermodynamics package [34]. In OPA, six primitive equations (momentum balance, the hydrostatic equilibrium, the incompressibility equation, the heat and salt conservation equations and an equation of state) are solved, where the Arakawa C grid is used in the horizontal. In the vertical, terrain-following coordinates, z coordinates or hybrid z-s coordinates can be chosen. Previously, NEMO was applied to the Baltic Sea and the North Sea area in uncoupled mode [35], coupled to atmospheric models [36] and forced with a wave model [18,19,25]. For the northwestern European Shelf, NEMO is used as a forecasting model in the COPERNICUS Marine Services [23,37,38]. The horizontal model resolution is about 2 nm with 51 σ-levels in the vertical providing instantaneous hourly surface (0.6 m uppermost level thickness) velocity fields. The study domain is 48.0–62.5° N and −4.7–13.2° E, which includes the North Sea, Skagerrak and Kattegat (Figure 1). The hourly atmospheric forcing is taken from subsequent short-range forecasts from the regional atmospheric model COSMO-EU, operated by the German Weather Service (DWD). Atmospheric pressure and tidal potential are included in the model forcing. River run-off is provided in the form of a daily climatology based on river discharge datasets. Lateral open boundary and initial condition fields (temperature, salinity, velocities and sea level) are derived from the MetOffice Forecasting Ocean Assimilation Model (FOAM) AMM7 (7 km horizontal resolution [23]), currently used by the Copernicus Marine Environment and Monitoring Service (CMEMS) as an operational service.

2.1.2. The Wave Model WAM

The wave model WAM [39,40] is a third-generation wave model, which solves the action balance equation without any a priori restriction on the evolution of spectrum. It is based on the spectral description of the wave conditions in frequency and directional space at each of the active model sea grid points of a certain model area. The version used in this study is the WAM Cycle 4.7, which is described in [41,42,43]. The source function integration scheme is made by [44], and the updated source terms of [45] are incorporated (Appendix A). The new version considers the wave-induced processes needed for coupling and are described in the next section. The spectrum in WAM is discretized with 24 directions and 25 frequencies. The wave model boundary information used at the open boundaries is taken from the regional wave model EWAM for Europe, which is run twice daily in operational routine at DWD.

2.1.3. Wave-Induced Processes

Ocean waves influence the circulation through a number of processes: turbulence due to breaking and non-breaking waves, momentum transfer from breaking waves to currents in deep and shallow water, wave interaction with planetary and local vorticity, Langmuir turbulence. The NEMO ocean model has been modified to take into account the following wave effects [11,18]: (1) the Stokes–Coriolis forcing; (2) sea-state dependent momentum flux, set as a scalar dependence of the flux from the atmosphere to waves and ocean or as a vector; and (3) a sea-state dependent energy flux. A description of the wave-induced forcing and the processes of wave interaction with the ocean circulation is given in Appendix A.1.

2.1.4. The Lagrangian Model

OpenDrift [46] is a freely available open-source off-line Lagrangian particle trajectory model that contains several modules for the advection of, e.g., oil spills [47], larvae and passive tracers [48]. For this study, the passive tracer module is used, which advects tracers only due to currents and winds. The sea and wind drift input is described in the next section. To investigate the influence of the wind, three different windage (referred to as leeway) coefficients when the bulk effect of waves and wind on the object is considered [32] (L_w) are used (0.1%, 0.5% and 1.0%). The real wind drag of the drifters is not known, so these values are used as estimates. The windage coefficients are multiplied by the wind velocity and added to the current velocity. The advection scheme is a 2nd-order Runge-Kutta scheme, and no additional diffusion is added because we only want to investigate the wave effects without any additional random “disturbances” [49]. Each particle thus follows a trajectory influenced by the sea surface currents and the equations read:

$$\mathsf{\Delta }x\left(t\right)=\left(u\left(t+\frac{\mathsf{\Delta }t}{2}\right)+L_w{u}_{wind}\left(t+\frac{\mathsf{\Delta }t}{2}\right)\right)\Delta t$$

(1)

$$\mathsf{\Delta }y\left(t\right)=\left(v\left(t+\frac{\mathsf{\Delta }t}{2}\right)+L_w{v}_{wind}\left(t+\frac{\mathsf{\Delta }t}{2}\right)\right)\mathsf{\Delta }t$$

(2)

here $\mathsf{\Delta }$x(t), $\mathsf{\Delta }$y(t) are the particle displacements and u(t + $\mathsf{\Delta }$t/2), v(t + $\mathsf{\Delta }$t/2) are the horizontal velocity components at the particle’s position at t + $\mathsf{\Delta }$t/2. If a particle reaches land, it is not further advected and considered beached. These particles are not further taken into account. If a particle leaves the domain through an open boundary, e.g., from the North Sea into the Atlantic, this particle is treated the same way. The model time step $\mathsf{\Delta }$t and the particle seeding strategies are dependent on the specific experiment (provided in Section 3).

2.2. Observational Data

2.2.1. Drifter Data

The HE 445 cruise was performed between May and July 2015 (see Figure 1b for the trajectories and deployment period of the drifters) [26]. The R/V Heincke deployed nine Albatros drifters corresponding to two models (Figure S1): MD03i (drifters 1–6) and ODi (drifters 7–9). The drifters provided their current positions by a Global Positioning System (GPS), which were transmitted to the R/V Heincke via Iridium (a bi-directional satellite communication network). The MD03i is a cylinder-shaped drifter, with a diameter of 0.1 m and length of 0.32 m, but only approximately 0.08 m remains above the water surface. The MD03 drag ratio is 33.2, and according to the parametrization [26] the MD03 slippage is around 1.1 to 1.6 cm/s, for 10 m/s wind speed and velocity difference across the drogue (ΔU) equal to 0.1 cm/s. The ODi is a spherical drifter with a diameter of 0.2 m, but only approximately 0.1 m is above the water surface. A sail (0.5 m in length and diameter, Figure S1) was attached to every drifter to enhance the drag below the water surface, and it was 0.5 m below the surface. Due to the small drifter surface above the water compared to the sail surface below the water, the drifter is designed to follow the ambient current in the upper meter of the water column. Due to the meteorological conditions, only some of the drifters were recovered at the end of the campaign. The recovered drifters corresponded to the short data-set, while the non-recovered drifters from the long data-set transmitted data until the batteries drained (Figure 1b).

2.2.2. HF Radar Data

HFR surface current data were acquired by three radar stations in the German Bight [50,51]. The radar systems are based on linear antenna arrays installed near the shoreline. More details on the required processing can be found in [52]. For the considered system, the spatial radar resolution is 1.5 km in range. The radar system in the German Bight can reach out to about 120 km off the coast in favorable conditions and provides measurements with a 9 min averaging window every 20 min. The observations are available as interpolated to a 2 km Cartesian grid. Through a combination of the radial components from the different antenna stations, meridional and zonal current components can be derived; however, the original radial components were used for the subsequent data assimilation procedure. Further details of the system can be found in [51].

2.3. Model Experiments

For the control experiment (REF), the velocity is taken from the stand-alone circulation model NEMO (Section 2.1), and the wave–current interaction processes are not included. In the coupled wave-ocean experiments (CPL), the wave-induced processes described in Section 2.1.3 calculated by the wave model WAM are introduced in NEMO to simulate the Eulerian velocity, needed for the OpenDrift. We performed the following experiments, in which the individual or combined effects of the wave-induced processes are included: (i) sea state dependent momentum flux (CPL-TAUOC); (ii) Stokes–Coriolis forcing (CPL-STCOR) and (iii) both the sea state dependent momentum flux and the Stokes–Coriolis (CPL-TAUST).

In order to study the role of the windage, two additional sets of experiments have been performed. For these, in addition to the Eulerian velocity from REF and CPL experiments, the windage is included in OpenDrift. We will name these experiments WD-REF and WD-CPL, respectively. The different windage contributions that we consider are 0.1, 0.5 and 1% (WD-REF and WD-CPL with _0.1, _0.5, and _1.0, respectively). All experiments are summarized in

Table 1. List of experiments.

Experiment	REF	CPL	WD-Ref	WD-CPL
NEMO-only	Yes	No	Yes
NEMO-WAM	No	Yes	No
Windage	No	No	0.1%/0.5%/1%	0.1%/0.5%/1%

.

3. Evaluation of Model Simulations

3.1. Methodology

The evaluation of the model results consisted of two parts. The first part (Section 3.2) aimed at studying the model runs statistically with direct comparisons of the simulations with the drifter data. In the follow-up part (Section 5), the drifter trajectories were modelled with OpenDrift to investigate the importance of the wave effects, e.g., for search-and-rescue applications.

The drifter velocities along the trajectories are calculated by dividing the distance between each drifter position by its time difference. The velocity is located in the middle between the two positions in time and space. Afterwards, the model velocities are interpolated trilinearly (lat, lon, time) to the drifter velocity positions and times. The root-mean-square error (rmse), standard deviation (std), bias (bia) and linear correlation coefficient (cor) are calculated to assess the model and drifter currents. Trilinear interpolation to the drifter position is performed for the wind and water depth and for each model experiment. For the analysis of meteorological conditions, periods with weak winds (25–27 June 2015) and strong southerly winds (01–03 June 2015) were chosen. In order to assign the errors of model velocities to different sources, a multi-linear regression was performed, solving the following equations:

$$\begin{array}{c}{u}_{drifter}=u_{model}{a}_{1i}+u_{wind}{b}_{1i}+c_{1i},\\ v_{drifter}=v_{model}{a}_{2i}+v_{wind}{b}_{2i}+c_{2i},\end{array}$$

(3)

Here, u_model is the model velocity, u_wind is the wind velocity, $U=\sqrt{u^2+v^2}$ denotes the velocity magnitude, i corresponds to the drifter number and a, b, c are the coefficients. Coefficient b represents direct windage acting on the GPS drifters. The coefficient a can indicate if the model under- or over-estimates surface velocity. Coefficient c includes deviations that cannot be explained by the model currents or the wind, like inaccurate directions; it can be considered an “offset”.

The drifter starting positions are used as initial positions in OpenDrift. In the trajectory simulations, the particles are initialized at the start location and start time of the Albatros drifters and move under the influence of the forcing until the real drifters expire. Drifter simulations of 25 h drift paths that were initialized every day from 0 to 53 at 13:00 UTC were previously performed [27]. It is important to stress that inall of our experiments, after starting the drifter simulations, we did not re-initialize the modelled drifters. The time steps in the OpenDrift experiments and the output frequency are both 10 min. In the first set of experiments, only the currents of REF and CPL were used in OpenDrift. To study the impact of the wind, REF and CPL current data were used together with the windage coefficients WD-REF and WD-CPL, respectively.

As a measure of the modelled trajectories, the skill score ss [4] is chosen. First, an index s is estimated as

$$s={\displaystyle \sum }_{i=1}^N{d}_i/{\displaystyle \sum }_{i=1}^N{l}_{oi}$$

(4)

where N is the total number of time steps, d_i is the distance of real and modelled drifter at time step i and l_oi is the total length of the trajectory of the real drifter at time step i. This index is used to calculate ss with

$$ss=\{\begin{array}{c}1-\frac{s}{n}, \left(s\le n\right)\\ 0, (s>n)\end{array}$$

(5)

where n is the tolerance threshold, equal to 1 [4], which means that the cumulative separation distance is not larger than the cumulative trajectory length. If s and ss are close to 1, the modelled trajectory is close to the observed trajectory. It assumes that the model performance is better with higher ss and the cumulative separation distance is less significant than the cumulative trajectory length. If the skill score is close to zero, the model simulations have no skill.

The skill of the particle transport model depends on the quality of the ocean circulation, but the skill score has limitations [53]; e.g., in situations with small currents, it yields too small cumulative distances. It can thus give a too high estimate of s and a low skill score ss. To overcome these limitations, it has been suggested [4] to have an appropriate definition of the tolerance threshold n. In the German Bight study area, due to the strong tidal currents and wind forcing, the skill score can be an applicable measure for the drift model performance.

3.2. Surface Currents

3.2.1. HF-Radar versus Drifter Observations

Figure 2 shows the drifter velocity magnitude versus the HF radar velocity as scatter plots generated for the whole period of drifter deployment. For the study period, 334 collocations of radar and drifter data were available for comparison. The velocity data of the drifters were filtered using a Savitzky–Golay filter [54] with a quadratic fit and a window length of ~2.3 days in order to remove outliers. The velocity of all nine drifters (regarding their deployment period or quality) was used in the scatter diagram to demonstrate the statistical robustness of the method. A detailed analysis of the separate trajectories of the Albatros drifters for the deployment periods was reported [27].

The HF radar currents should be interpreted as Eulerian currents [55] (i.e., without including the Stokes drift). However, it has also been argued that the HF radar velocity partially contains the Stokes drift, in which case the HF radar measurements contain a filtered component of the Stokes drift [56]. On the other side, comparisons of HF radar velocity measurements with drifter observations demonstrated that the presence of the Stokes drift in HF measurements is not settled. In Figure 2 and Figure S2 we compare HF radar velocities interpreted as Eulerian currents against the velocity estimated by the nine Albatros drifters. With a standard deviation of 15 cm/s and a bias of 2.7 cm/s, the drifter and the HF radar are in reasonably good agreement. A noticeable feature is a slight underestimation of HF velocity at higher values. There are several reasons for the disagreement between HF and drifters when the speed is over 0.6 m/s: (i) the HF radar retrieves the currents at a depth of around 1 m, but the drifter’s sails are located at 0.5 m. An underestimation of the HF radar regarding the drifter measurements may be expected. (ii) Due to strong wind stress during storm events, taking into account that the current profile over depth is logarithmic, a stronger current gradient between 1 m (HF) and 0.5 m (drifters) depth is expected. (iii) During breaking events in storm conditions, the drifter and the sail can surf over the sea surface, while the HF radar measures the current speed under the surface.

3.2.2. Assessment of Model Velocity

The simulated model velocities of the stand-alone NEMO (REF) and wave-circulation coupled model (CPL) are compared with the matching drifter velocity in Figure 3a. REF and CPL model velocities are well represented in the scatter plots and show good agreement with the drifter data. Above 30 cm/s, the modelled velocities are lower than the drifter velocities, as seen in the zonal and meridional components. The range of the q-q plot is smoother and fits better in moderate and high velocities in the coupled model simulations (Figure 3). The comparisons demonstrate a better fit of the meridional velocity to the observations by the coupled wave–current model for strong winds in both directions (Figure S3). The Taylor diagram shows a good agreement between both simulations and observations with slightly improved skills of the CPL experiments. It is noteworthy that despite the coarse resolution of the model and the complex bathymetry and coastline in our study area, the comparison is satisfying, and a general improvement of the surface currents in both directions is observed due to the wind–wave–ocean coupling. This is consistent with previous findings [21] in which coupled and stand-alone model circulation against Acoustic Doppler Current Profiler (ADCP) demonstrated an intensification of velocities due to coupling with waves, leading to a reduction in prediction errors.

A multi-linear regression is performed, including model velocities, wind velocities and an offset to obtain insight into the error sources. Figure 4 shows these coefficients separately for each of the nine drifters. The model velocities are reduced in the zonal component by about 2.5% and in the meridional component by about 8%. While reducing the model velocities, the wind to be taken into account is around 1.0% for both velocity components. This value can also be interpreted as a guess about the wind drag of the drifters. For the same drifters, different windage and Stokes drift contributions were tested [27], in a combination of the Eulerian velocity of two hydrodynamic models. The parameterization used in [57] predicted a slippage of 1.1 to 1.6 cm/s for 10 m/s winds, and the expected windage should thus be close to 0.135%. However, as demonstrated in [58] this value depends strongly on the drifters and the chosen models and usually ranges between 0.1 and 1%. A wind drag of 0.27% was calculated from drifter surface ratio and 0.3% with Stokes drift was found to be the best combination in model simulations [49]. The offset of the zonal components is slightly negative (about −2 cm/s), while the offset of the meridional ones is slightly positive (about 2 cm/s).

By considering the wave coupling in the circulation model (Figure 4, right panels), the windage halves and model velocities are taken less into account. Instead, the offset increases. Note that drifter 7 shows substantial differences compared to the others, which may be due to technical problems. Drifters 1 and 9 also show increased/decreased coefficients. Together with drifter 7, these are the only drifters that are beached. A possible explanation is that the drogue was damaged during the storm and then landed on the beach. Thus, model uncertainties at the boundaries or grounding of drifters could be a possible reason.

Figure 5 shows the time series for the chosen periods of drifter velocity and wind conditions for the stand-alone NEMO model and coupled model experiments considering wave-induced processes described in Section 2. The tidal cycle is simulated relatively well by all model experiments. There is a good fit between the observations and the model velocities during calm wind and wave conditions. The magnitude of the velocity in the CPL-TAUOC experiments is higher on 1 June than in REF, due to the veering of the wind, which influences the sea-state dependent momentum flux forcing. In periods of strong winds (e.g., 03/06), the velocity of the REF model is under-estimated, while the coupled NEMO-WAM currents performed better than the currents from the stand-alone NEMO mode. The wave model compares very well against in situ and satellite observations [43,59], see Figure S4.

4. Model Trajectories

We performed several Lagrangian experiments to investigate the impact of wave-induced forcing on the ocean model. In the first set of experiments, OpenDrift used only the Eulerian current fields as provided by NEMO (REF and CPL experiments) without adding any wind drift correction. In the second set of experiments, an additional windage was included.

4.1. Time Series of Separation Distance, Skill Scores and Standard Deviation

As a separation metric, the normalized cumulative Lagrangian separation [4] was applied. First, we focused on the model skill in drifter trajectories to demonstrate the impact of the local wind and wave-induced velocity correction terms. Further, we quantified the sensitivity of the simulated particles to the individual or combined effect of the wave-induced process that were implemented in our coupled wave-circulation model system.

During the whole integration period, the separation distance of the drifters (Figure 6) remained very low, even though we did not perform any re-initialization of the drifter model after the initial release. Drifters 2–4 travelled for only a few days in the German Bight (Figure 1). It is noticeable that during the lifetime of these drifters, the simulated separation distance was kept within the model grid resolution. The model performance for drifters 5 and 6 (Figure S5) was high for all coupled experiments, and the separation distance was kept below 20 and 40 km, respectively, even one month after the start of integration. The separating distance of the REF (drifters 5–6) sharply increased after 20 June, showing significant deviation from the observations. This result coincides with earlier findings [27], demonstrating that on 16 and 22 June the wind direction sharply changed. Consequently, under these transitional conditions, signs of directional errors substantially differed, and the model performance was unsatisfactory. Implementing wave-induced processes into NEMO and using these Eulerian currents in OpenDrift led to a decrease in the separation distance between the CPL model experiments and the observations. In [18] was found that intensification of zonal velocity towards the coast, simulated by the coupled model, fits better with ADCP measurements. The present results show that for the CPL experiments, even during periods with moderate significant wave height, the inclusion of wave parameterizations improved the model performance. This demonstrates that separately adding a contribution of the Stokes drift to the Eulerian currents to OpenDrift in many cases might be insufficient to simulate the drifters under changing sea state conditions appropriately. The external Stokes drift can also be inconsistent with simulated ocean currents. For drifter 3 (Figure 6a), the distance error was of the same order as the resolution of our ocean model. For integration periods longer than a month, the distance error was less than twice the order of the grid resolution for the one-month simulation of drifter 5 (Figure 6c). Thus, the sea-state contribution of the momentum flux helps bring down the separation distance during most of the simulation periods.

Sensitivity experiments on different wind drift (leeway) percentages (Figure 7 and Figure 8) demonstrated that only very few of them managed to reach higher skill scores than the CPL experiments. The results show the significance of producing the Eulerian velocity by the circulation models forced by sea state dependent momentum fluxes. The skill scores of the simulated drifters were generally high in all experiments, above 0.8 (Figure 9). The CPL experiments managed to keep the std low even after this period without considering any windage (see Figures S6–S8).

Applying wind drift correction of 0.5% or higher to WD-CPL experiments reduced the skill score (Figure 9c,d). In the first days of the integration period, the skill score of drifter #5 by windage correction of 0.5% and 1.0% dropped to 0.7 and 0.5, respectively. The skill score of the experiment with windage correction of 0.1% was closer to the CPL experiments but still slightly lower in the period from 26 May until 2 June (see also Figures S7 and S8). These differences were illustrated even better by the standard deviations of the WD-CPL experiments (Figure 9d). The standard deviations of D#5 stayed almost constant through the integration, at about 0.5 and 5 km, respectively, for CPL and WD-CPL plus 0.1% (Figure 9). Adding 0.5% windage yielded a higher standard deviation for all experiments. By assuming 1.0% windage, a trend of the systematic increase was observed for all drifters. Neither of the WD-CPL experiments (in Figure 9c,d) improved the skill scores and reduced the deviation between the trajectories of the Lagrangian model and observations by additionally considering the windage contribution to the Eulerian currents obtained by the coupled wave–circulation ocean model.

Sensitivity to the windage contribution was also determined for the REF runs (Figure 9e,f). The skill score and standard deviation of WD-REF were better than those of WD-CPL for all experiments. On the other side, the skill scores/standard deviations were lower/higher than those of CPL. Only by considering 0.5% windage of WD-REF are the ss and dd similar to the CPL.

It was demonstrated in [53] that simulated trajectories obtained by considering the wind drag of current velocities or the Stokes drift’s contribution showed better skills against the observations than without the corrections. They postulated that considering 1.0% windage or adding the contribution of the Stokes drift gives almost identical results. Our sensitivity analysis, however, showed that the CPL experiments performed best. From the rest of the experiments, only WD-REF with 0.5% windage provided similar skill. The model experiments showed that using Eulerian velocity estimated by considering Stokes–Coriolis and the sea state momentum provided similar results.

4.2. Particle Trajectories of the Albatros Drifters

The trajectory of drifter 5 in the REF experiment deviated in a northwesterly direction. (Figure 10a) The trajectories obtained by the CPL experiments correctly reproduced the drifter direction and were in good agreement with observations. By north and northwesterly winds on 10 June and the anticyclonic circulation [27], the CPL experiments reproduced the windage of drifter 5 well. CPL-TAUOC separation distance was on the order of the model grid resolution (Figure 6). The low dd values remained until 29 June. The trajectory of WD-CPL_1.0 (Figure 10b) made a higher loop towards the north and northeast, deviating from the observations. By a 0.5% contribution of windage, the drifter was unrealistically transported towards the east and was beached. The trajectory of WD-CPL with 0.1% started deviating for 2.3 days after starting the simulation. WD-REF with 0.5% leeway showed a smaller separation distance (Figure 8), and the trajectory remained closer than that of WD-CPL with 0.5%. By adding 1.0% leeway, the simulated trajectories of WD-CPL and WD-REF can be considered as wrong. Skill scores and STDs for all experiments and during the deployment period of the drifters are given in Tables S1 and S2.

5. Discussion

Coupled ocean–wave models together with a Lagrangian transport model can be beneficial for drift and transport studies, ranging from search and rescue modelling of specific objects [33], backtracking [60], larval drift [61] and marine plastics [62,63] to the connectivity between different marine protected areas [16]. Accurate measures can potentially have a strong impact on biodiversity risk assessments (e.g., connected to marine litter or oil spills). In [16], several models of the North Sea were compared to study the variability and resulting uncertainty and differences between the models.

In our study, Eulerian currents from a coupled ocean-wave model were used to perform CPL and WD-CPL experiments. In this way, additional tuning of the contributing factors of the Stokes drift by the drifter model can be excluded in CPL experiments. Our results are in line with the study of [56]. The latter concluded that implementing the Stokes drift as a simple additive component of drift velocity, parameterized in terms of wind forcing, can be inconsistent (a violation of momentum and energy conservation) if Eulerian currents were simulated without taking into account the reservoir of wave momentum and energy. Our results indicate the relevance of the role of waves for redistribution of momentum, especially in periods of changing wind direction (Figure 6, Figure 7 and Figure 8).

The Eulerian currents in most of the Lagrangian transport models, e.g., for search and rescue operations, are taken from, for example, operational model output. The wind stress is parameterized to drive the dynamics of the upper ocean directly if the wave model is not included. However, part of the wind stress is supported by the flux of momentum from wind to waves. These processes were considered in our CPL-TAUOC and CPL-TAUST experiments to simulate the Eulerian velocity needed for OpenDrift. Due to the non-linearity of wave–current interactions, the individual effect of the wave-induced coupled STCOR and TAUOC processes are not superimposed on TAUST [18]. This result shows that the Eulerian currents by the coupled run provided the best fit to the observed particle trajectories. The worsening of the Lagrangian model skill, especially at the beginning of the drifter simulation, can have an impact on applications like search and rescue, in which skilled Lagrangian transport forecasts are needed at the very beginning of the operation. In summary, the use of additional wind corrections cannot be justified in WD-CPL experiments, since neither of the additional windage correction experiments demonstrated an improvement compared with the CPL experiments.

Displacements of the more offshore drifters, 5 and 6, were observed on 3–6 June [27], and the models turned out to be largely overestimated, concluding that neither simulated currents, windage fields nor Stokes drift by their Lagrangian model were able to reproduce the spatial gradients. In the present study, CPL experiments showed good agreement for these drifters with the observations. A deviation of the trajectory, and consequently an increase in the separating distance, was observed, but stayed within the same grid resolution. On the other side, the WD experiments deviations from the simulated drifter trajectories from the observations and separating distances were high. This result proves the importance of improving the Eulerian current simulations needed for Lagrangian transport modelling by implementing the sea state momentum forcing and Stokes–Coriolis forcing in the numerical model.

We did not aim to assess all specific differences between individual drifters and discuss in detail the wind forcing and the North Sea circulation during the different periods of their deployment, as this was done by [27]. Neither did we aim to tune the model to find the optimal percentage of the contribution of either the Stokes drift or windage that can be taken into account in OpenDrift simulations (as was recently done by [64]. Additional errors of the drifter simulations can be due to the errors in the atmospheric forcing, but also the boundary or tidal forcing. Another reason may be errors in the vertical resolution and mixing and bottom friction parameterization; the latter is essential for shallow water dynamics in regions like the German Bight. By using the same model configuration as in the present work, the interactions between barotropic tides and mesoscale processes were studied by [65] showing that barotropic tides affect diapycnal mixing substantially. In our study, we aimed to investigate the potential of a different approach, namely the contribution of wave-induced processes in the NEMO model.

It is important to stress here again that the North Sea circulation is very complex. The drifters were deployed in the shallow German Bight, a coastal ocean where the currents are dominated by tidal and wind forcing and are steered by the bathymetry and coastline. Other possible error sources are circulation features such as inertial oscillations, sub-mesoscale dynamics and baroclinic effects. Further studies are needed to quantify the combined effects of the mesoscale variability and resolution of the ocean and particle transport model, especially in the coastal areas of the German Bight (the Wadden Sea) as well as for other regions. For these aims, the GCOAST framework will be taken into account and tested for its best performance in terms of trajectory simulations over regions with different oceanographic properties (e.g., Baltic Sea, the North-East Atlantic).

The advantage of the method here is that by coupling the circulation model to the wave model, no further sensitivity experiments on the percentage of the contribution of windage or Stokes drift to improve the skills of the Lagrangian model are needed.

6. Conclusions

The results of our experiments lead us to the following conclusions:

1. Comparing currents from coupled and stand-alone model simulations demonstrated that the coupled model velocities fit better with the observations, especially moderate and high values. In calm wind and wave conditions, the differences are not pronounced. By using a fully coupled model, consistent atmosphere–wave–ocean forcing is applied to simulate the Eulerian currents needed for particle transport simulation. Besides, the bias of the directions in the wind–wave–ocean currents simulations is reduced compared to that in the stand-alone model. It is noteworthy that despite the coarse resolution of the NEMO model and the coupled NEMO-WAM setup and the complex bathymetry and coastline in our study area, the comparisons are satisfactory, also for the stand-alone model simulations. A general improvement in surface currents for both directions is observed due to the wind–wave–ocean coupling.

2. The multi-linear regression analysis showed that for CPL, windage is halved (from 1.0 to 0.5%), and the model velocities are taken less into account. Instead, the offset increases. These results show that the wind drift is better accounted for in the coupled NEMO-WAM model than in the stand-alone NEMO. It was also shown that the regression can reveal technical problems of the drifter, like probable drogue loss.

3. We showed the particle analysis by calculating values such as separating distance, skill score and standard deviation between model experiments and observations. We demonstrated that the skill score, based on the cumulative Lagrangian separation distance standardized by the associated cumulative trajectory length, proved to be a useful parameter to evaluate the overall model performance, rather than using a daily validation metric. Although the skill score and standard deviation of the ODi drifter are slightly lower than those of the MD03i, both types demonstrate good predictive skill. The MD03i drifter shows better skill for CPL than ODi. The skill scores with wind showed similar behavior by both types of drifters. For drifters that reached land, the model performance was low, which might be due to the insufficient model resolution and parameterization together with higher model uncertainty at the boundaries or beaching of drifters.

4. The skill scores and separation distances improved when wave-induced processes were taken into account in the ocean-only simulations. By considering sea state momentum dependencies or Stokes–Coriolis forcing in the coupled model, the skill scores are quite similar. Adding the contribution of windage or Stokes drift to currents produced by the fully coupled run with waves might lead to over-parameterization.

5. In the CPL experiments, it turned out that no additional tuning of the wind drift factor was needed for best fit with the observations. This also indicates that no additional contribution of windage or separate consideration of external Stokes drift was needed to predict the drifter trajectories satisfactorily. For some drifters, the skill scores of the WD-REF or WD-CPL experiments were similar to the CPL by adding a direct windage contribution to the REF. However, these percentages varied between different drifters and experiments.

The results based on the drifters used in this study showed that, in some cases, it might be favorable to implement full coupling of waves and circulation models to produce the currents needed for drifter simulations. This could lead to a more consistent approach instead of trying to tune the windage factor or percentage of external Stokes drift, separately or combined. Such tuning of the contribution coefficients is typically restricted by the need for availability and testing a large number of drifter observations (that are normally lacking), which is required to improve wind–wave–current forcing dependencies in the Lagrangian model. We note, however, that this conclusion is based on a limited period and a small area (over the German Bight). Nevertheless, the model simulations showed that the results are promising for better understanding and prediction of Lagrangian transport by using Eulerian currents simulated by coupled wave–ocean model simulations instead of a stand-alone ocean model. The results revealed that the newly introduced wave effects are essential for the drift model performance.