Empirical Lagrangian parametrization for wind-driven mixing of buoyant particles at the ocean surface

Turbulent mixing is a vital component of vertical particulate transport, but ocean global circulation models (OGCMs) generally have low resolution representations of near-surface mixing. Furthermore, turbulence data is often not provided in OGCM model output. We present 1D parametrizations of wind-driven turbulent mixing in the ocean surface mixed layer, which are designed to be easily included in 3D Lagrangian model experiments. Stochastic transport is computed by Markov-0 5 or Markov-1 models, and we discuss the advantages/disadvantages of two vertical profiles for the vertical diffusion coefficient Kz . All vertical diffusion profiles and stochastic transport models lead to stable concentration profiles for buoyant particles, which for particles with rise velocities of 0.03 and 0.003 m s−1 agree relatively well with concentration profiles from field measurements of microplastics when Langmuir-circulation-driven turbulence is accounted for. Markov-0 models provide good model performance for integration timesteps of ∆t≈ 30 seconds, and can be readily applied in studying the behaviour of 10 buoyant particulates in the ocean. Markov-1 models do not consistently improve model performance relative to Markov-0 models, and require an additional parameter that is poorly constrained. Copyright statement.

as: (1) where K z = K z Z(t) is the vertical diffusion coefficient, ∂ z K z = ∂K z /∂z, dW is a Wiener increment with zero mean and variance dt and we define the vertical axis z as positive upward with z = 0 at the air-sea interface. The Euler-Maruyama (EM) 60 scheme (Maruyama, 1955) is the simplest numerical approximation of equation 1, where infinitesimal terms dt and dW are replaced with the finite ∆t and ∆W . Equation 1 can then be rewritten as (Gräwe et al., 2012): where w is the stochastic velocity perturbation due to turbulence. The turbulent transport has both a deterministic drift term the autocorrelation of the particle velocity from its initial velocity due to unresolved sub-grid processes, which depends on the model resolution and setup in a given study. Since there is not a clear indication of the true value of T L , we consider a range of values α ∈ [0, 0.1, 0.3, 0.5, 0.7, 0.95], corresponding to T L ∈ [1, 1.1, 1.4, 2, 3.3, 20] × ∆t. As the depth dependence of T L is uncertain, we make the simplification that ∂ z T L = ∂ z α = 0. Since ∆t ≤ T L , we use K z = σ 2 w ∆t (Brickman and Smith, 2002), 90 which means that equation 6 becomes: In this form, it is clear that equation 7 is equivalent to equation 4 when α = 0. This is because when α = 0, velocity perturbations w are assumed to be uncorrelated over timescales ≥ ∆t, which is equivalent to the M-0 formulation. M-1 stochastic models generally should lead to improved representation of diffusion in Lagrangian models (Berloff and McWilliams, 2003;95 Van Sebille et al., 2018), but it does require insight into turbulence statistics that have not yet been extensively studied in Lagrangian settings. For that reason, while even higher order Markov models are theoretically possible (Berloff and McWilliams, 2.2 Vertical diffusion profiles 120 Two vertical diffusion coefficient profiles are used, with the first based on Kukulka et al. (2012) and Poulain (2020). Kukulka et al. (2012) parametrized the near-surface vertical diffusion coefficient K S z due to breaking waves as: for z > −1.5H s , where κ = 0.4 is the von Karman constant, H s is the significant wave height and u * w is the friction velocity of water. The significant wave height H s is parametrized as H s = 0.96g −1 β 3/2 * u 2 * a , where g = 9.81 m s −2 is the accelation 125 of gravity, β * = c p /u * a is the wave age, c p being the characteristic phase speed of the surface waves and u * a = τ /ρ a is the friction velocity of water. The friction velocity of air is based on the air density ρ a = 1.22 kg m −3 and the surface wind stress τ = C D ρ a u 2 10 , where u 10 is the 10m wind speed and C D is the drag coefficient (Large and Pond, 1981). Similarly, u * w = τ /ρ w with the seawater density ρ w = 1027 kg m −3 . Following Kukulka et al. (2012), we assume a fully developed sea-state with β * = 35. The Kukulka et al. (2012) parametrization is valid only for z ≈ −1.5H s , and we extend the parametrization for greater 130 depths using the eddy viscosity profile ν z as found for oscillating grid turbulence by Poulain (2020): where ν S is the near surface eddy viscosity. This approach agrees with Kukulka et al. (2012) in predicting constant mixing for z > −H s , where the eddy viscosity then drops proportional to z −3/2 for greater depths. Oscillating grid turbulence (OGT) experiments are commonly used to study wave and wind induced turbulence (Fernando, 1991). As OGT experiments have 135 been shown to reproduce turbulence decay laws of velocities and dissipation rates observed in the ocean ML (Thompson and Turner, 1975;Hopfinger and Toly, 1976;Craig and Banner, 1994), this provides some confidence in the modelling of the decay of near-surface eddy viscosity, although direct validation with field measurements of eddy viscosity have yet to occur. The diffusion coefficient K z depends on ν z as K z = ν z /Sc t , where Sc t is the turbulent Schmidt number, and assuming ∂ z Sc t = 0, combining equations 8 and 9 results in: where K B = 3 × 10 −5 m 2 s −1 is the dianeutral diffusion below the MLD (Waterhouse et al., 2014). The diffusion is thus constant for z > −H s , below which K z ∝ |z| −3/2 , while the magnitude of K z increases for higher wind speeds (Fig. 1). As z → −∞, |z| −3/2 → 0, and therefore we include the bulk dianeutral diffusion K B to account for vertical mixing at depths below the influence of surface wave-driven turbulence. As both Kukulka et al. (2012) and Poulain et al. (2018) considered tur-145 bulence generated by breaking surface waves, we refer to this diffusion approach as Surface Wave Breaking (SWB) diffusion.
The second vertical diffusion coefficient profile is a local form of the K-profile parameterization (KPP) (Large et al., 1994; where φ = 0.9 is the "stability function" of the Monin-Obukov boundary layer theory, θ is a Langmuir circulation (LC) enhancement factor, and z 0 is the roughness scale of turbulence. As such, K Z rises from a small non-zero value at z = 0 to a maxima at z = 1/3M LD, before dropping to K z = K B for z ≤ M LD (Fig. 1). In the original KPP formulation K z (z ≤ M LD) = 0 since the theory only applies to the surface mixed layer, so we add the same bulk dianeutral diffusion term K B as with the and has been shown to strongly affect the vertical concentration profiles of buoyant microplastic particles in LES experiments (Brunner et al., 2015;Kukulka and Brunner, 2015). Therefore, we examine θ ∈ [1.0, 2.0, 3.0, 4.0, 5.0]. The roughness scale z 0 , which can represent the surface roughness due to surface waves, depends on the wind speed and the wave age (Zhao and Li, 2019), and following Kukulka et al. (2012) we consider a wave age β * = c p /u * a = 35 that is equivalent to β = c p /u 10 = 1.21.

160
According to Zhao and Li (2019), the roughness scale is given by: For w 10 = 0.85 − 9.30 m s −1 , this means z 0 = 2.38 × 10 −6 − 2.86 × 10 −4 m. To test the model sensitivity to z 0 , we also consider an alternative scenario where z 0 = 0.1 × H s = 1.76 × 10 −3 − 2.10 × 10 −1 m, following the same formulation H s = 0.96g −1 β 3/2 * u 2 * a as in Kukulka et al. (2012). This increases K z for z ≈ 0, but does not significantly affect the magnitude K z 165 at greater depths ( Figure B1). The original KPP theory does not explicitly account for surface wave breaking, which would lead to larger non-zero K z at z = 0. While we do not claim that setting z 0 = 0.1 × H s means that our KPP profile accounts for surface wave breaking turbulent mixing, it allows us to investigate the influence higher near-surface mixing would have on the modelled vertical concentration profiles. The MLD is the maximum depth of the surface ocean boundary layer formed due to interaction with the atmosphere, and in KPP theory the MLD is defined as the depth where the bulk Richardson number Ri B

170
is first equal to a critical value Ri crit . In the original formulation Ri crit = 0.3 (Large et al., 1994), but Ri B can be difficult to compute in the field as this requires data for both vertical density and velocity shear profiles. In this study we prescribe MLD= 20 m, as this falls within the range of the MLD for field data used to evaluate the model (see Section 2.3).

Field data
We compiled a dataset of vertical plastic concentration profiles collected within the surface mixing layer to validate the mod-175 elled concentration profiles (Table 1)

185
Almost all measurements were collected with neuston nets, either multi-level nets simultaneously sampling fixed depth intervals (Kooi et al., 2016b) or using multi-stage nets that consecutively sample fixed depths or depth ranges (Kukulka et al.  All measured microplastic concentrations are normalized by total amount of plastic measured within a vertical profile. In 195 order to compare the average normalized field concentration with the modelled profiles, we bin the normalized field concentrations into 0.5 m depth bins and calculate the standard deviation for each depth bin. Comparison of the modelled concentration profiles with the binned normalized field measurements is done via the root mean square error (RMSE): where C f,i and C m,i are the binned normalized field measurement and modelled concentration within depth bin i. Model 200 evaluation for the low buoyancy particles is not possible with the available field measurements as low buoyancy particles are typically too small to be sampled with neuston nets, and the Pieper et al. (2019) dataset alone is too small.

Results
Starting with all particles at z = 0 for t = 0, M-0 models with both KPP and SWB diffusion lead to stable vertical concentration profiles (Fig. 2), where the equilibrium concentration profile is already established within 1 -2 hours (Fig. C1). For 205 both diffusion profiles, there is progressively deeper mixing of particles with increasing wind speeds and decreasing buoyancy.
While with both SWB and KPP diffusion low buoyancy particles always get mixed below the surface, for medium and high buoyancy particles there exist minimum wind speeds below which all particles remain at the surface. These limits are similar for both diffusion types for medium buoyancy particles (u 10 ≥ 2.40 m s −1 ), but high buoyancy particles only mix below the surface with SWB diffusion if u 10 ≥ 9.30 m s −1 . However, once mixing below the ocean surface occurs, KPP diffusion always 210 leads to deeper mixing of particles than SWB diffusion due to higher subsurface K z values.
The concentration profiles for medium and low buoyancy particles are largely unaffected by reducing ∆t below 30 seconds ( Fig. E1). However, for high buoyancy particles with SWB diffusion the concentration profile more strongly depends on ∆t due to the applied boundary condition. For ∆t = 30 s, the M-0 model shows all particles remain near the ocean surface, but shorter 215 ∆t values indicate that deeper mixing of particles already occurs for u 10 = 6.65 m s −1 . With KPP diffusion, all high buoyancy particles remain at the surface even with ∆t = 1 second, as K z at z = 0 remains too low to overcome the high rise velocity.
Even though KPP diffusion with θ = 1.0 and z 0 following (Zhao and Li, 2019) predicts deeper mixing of particles than with SWB diffusion, both approaches underpredict the mixing of particles relative to field observations. For KPP diffusion, this can 220 be corrected by accounting for LC-driven mixing, which leads to deeper mixing of particles for both medium and low buoyancy particles (Figures 3 & D1). For medium buoyancy particles this generally leads to better model agreement with lower RMSE values between the modelled and averaged field data concentration profiles (Figure 4). However, for high buoyancy particles LC-driven circulation is not enough as particles remain at the ocean surface for all wind conditions even for θ = 5.0 ( Figure   D2), as K z for z ≈ 0 is too low to overcome the inherent particle buoyancy. Only when LC-driven is combined with higher 225 near-surface K z values by setting z 0 = 0.1 × H s do we see any below-surface mixing of high buoyancy particles when θ > 3.0 and u 10 ≥ 9.30 m s −1 . Increased near-surface K z values have a lesser influence on the concentration profiles of medium and low density particles, as these particles were already being mixed below the surface even without larger z 0 values.
With both KPP and SWB diffusion, M-1 models show deeper mixing of particles as α → 1 (Fig. 5). Relative to the field measurements, M-1 models can at best slightly improve model performance over M-0 models (Fig. 6). However, improved 230 model performance is not shown across all particle sizes and wind conditions, and there is not a consistent α value leading to the smallest RMSE values.

Discussion
The parametrizations presented in this study are intended for use in 3D Lagrangian experiments using OGCM data, and therefore should yield numerically stable results for the relatively large integration timesteps used in large-scale Lagrangian vertical 235 transport modelling (Lobelle et al., 2021). While there are more stable schemes available than the EM scheme used in this study (Gräwe et al., 2012), the EM scheme is computationally the cheapest and yields concentration profiles that match reasonably well with observations. Both M-0 and M-1 models show largely convergent concentration profiles for ∆t = 30 seconds, which would make both approaches feasible with regards to computational cost. However, we would currently recommend using a M-  acceptable. Otherwise, the error can be reduced by using a smaller integration timestep where that is computationally feasible.
Considering the KPP and SWB diffusion profiles, the results in this study indicate that KPP diffusion generally performs better relative to field observations. For high buoyancy particles, SWB diffusion leads to slightly deeper particle mixing, while only if the KPP diffusion profile accounts for LC-driven turbulence and has higher near-surface K z values can it similarly show 265 below-surface mixing of high buoyancy particles for u 10 ≥ 9.30 m s −1 . However, with medium and low buoyancy particles the KPP profile leads to much deeper mixing, especially when accounting for LC-driven turbulence, and this appears to agree better with field observations. As with (Brunner et al., 2015), we see that elevated near-surface K z values due to e.g. wave breaking have a comparably smaller influence on the overall concentration profile than LC-driven mixing, as similarly shown by (Brunner et al., 2015). Therefore, although we recommend future work incorporating surface wave breaking into KPP the-270 ory, our current KPP diffusion approach representing LC-driving mixing through θ seems to capture the majority of turbulent mixing dynamics. It must be noted though that the majority of the field measurements are collected in the top 5 meters of the water column, and more measurements would need to be collected at greater depths to further evaluate the depth to which  buoyant particles in LES model studies (Liang et al., 2012;Yang et al., 2014;Brunner et al., 2015;Taylor, 2018). The modelled concentration profiles generally resembled the profiles from field measurements of microplastic concentrations under different wind conditions (Kooi et al., 2016b;Kukulka et al., 2012), but the averaged concentration profiles of the field measurements are quite noisy. Partly, this could be due to inhomogeneity in the particle buoyancy, as the collected microplastic particulates 295 have varying sizes and rise velocities (Kooi et al., 2016b;Egger et al., 2020). Additionally, we sorted the field measurements based on wind conditions, but other underlying oceanographic conditions such as the MLD can still vary significantly even with similar wind speeds. Unfortunately, we lack additional data of the oceanographic conditions at the of sampling, which currently prohibits more high-level comparisons of the field and model concentration profiles. Compared with the field data, the variance in the modelled concentration profiles is significantly smaller. This is in part also due to assuming constant environ-300 mental conditions over 12 hours for the model simulations, while e.g. wind conditions can change on much shorter timescales over the ocean surface. To further improve vertical transport model verification, more measurements would be required, covering a wider range of oceanographic conditions (such as for wind conditions higher than u 10 = 10.7 m s −1 ) and with a high spatial sampling resolution also for depths z < −5m. Ideally these measurements would also sample small, neutrally buoyant particulates, but we acknowledge this is difficult with the sampling techniques commonly used today. At the same time, we 305 would encourage conducting more ocean field measurements of near-surface vertical eddy diffusion coefficient and/or eddy viscosity profiles, as this will allow further validation of the K z profiles predicted by the KPP and SWB theory with actual ocean near-surface mixing measurements.
The parameterizations have been validated for high/medium rise velocities, and at least for KPP diffusion with θ > 1.0 the 310 concentration profiles resemble those calculated from field observations. This provides confidence in the turbulence estimates from the KPP approach, and as these are independent of the type of particle that might be present, this would suggest the KPP approach can also be applied to to neutral or negatively buoyant particles. However, as model verification was only possible for microplastic particulates with rise velocities approximately between 0.03 -0.003 m s −1 , we would advise additional model verification for other particle types where the necessary field data is available. In the case of SWB diffusion, turbulent mix-315 ing seems underestimated when further from the ocean surface, and we would advise more validation with field observations before applying this diffusion approach to other particle types.
We have developed a number of 1D surface-mixing parametrizations designed to be readily applied in large-scale oceanic Lagrangian model experiments using OGCM data. Where possible, we would recommend using the turbulence fields from the 320 OGCM to assure turbulent transport of the particles is consistent with that of other model tracers. However, if the turbulence fields are unavailable then particularly parametrizations with KPP diffusion with LC-driven mixing are shown to produce modelled vertical concentration profiles that match relatively well with field observations of microplastics. The parametrizations generally perform well for timesteps of ∆t = 30 seconds, but for high buoyancy particles users need to take care to use sufficiently short timesteps, especially with SWB diffusion. Verification was only possible for positively buoyant particles larger 325 than 0.33 mm (which generally have rise velocities ≤ 0.003 m s −1 ), but the parametrizations should also be applicable to other particle types. The parametrizations can therefore be applied to investigate the influence of turbulent mixing on the vertical transport of (microplastic) particles within a 3D model setup, and ultimately gain a more complete understanding of the fate of such particles in the ocean.
6 Code and data availability