Integrating CVMix into GOTM (v6.0): A consistent framework for testing, comparing, and applying ocean mixing schemes

The General Ocean Turbulence Model (GOTM) is a one-dimensional water column model including a set of stateof-the-art turbulence closure models, and has widely been used in various applications in the ocean modeling community. Here we extend GOTM to include a set of newly developed ocean surface vertical mixing parameterizations of Langmuir turbulence via coupling with the Community Vertical Mixing Project (CVMix). A Stokes drift module is also implemented in GOTM to provide the necessary ocean surface waves information to the Langmuir turbulence parameterizations, as well 5 as to facilitate future development and evaluation of new Langmuir turbulence parameterizations. In addition, a streamlined workflow with Python and Jupyter Notebook is also described, enabled by the newly developed and more flexible configuration capability of GOTM. The newly implemented Langmuir turbulence parameterizations are evaluated against theoretical scalings and available observations in four test cases, including an idealized wind-driven entrainment case and three realistic cases at ocean station Papa, the northern North Sea and the central Gotland Sea, and compared with the existing General Length Scale 10 scheme in GOTM. The results are consistent with previous studies. This development extends the capability of GOTM towards including the effects of ocean surface waves and provides useful toolsets for the ocean modeling community to further study the effects of Langmuir turbulence in a broader scope. Copyright statement. TEXT

will allow us to precisely separate differences induced by variations in the turbulence parameterizations from those associated with different numerical schemes, parameterizations of the ocean surface fluxes, and other secondary effects.
In addition to documenting the implementation of CVMix in GOTM and the updates to improve the user interface, this paper also evaluates these new OSBL parameterizations in GOTM against available observations at a few sites. In particular, we focus on the influence of Langmuir turbulence in such evaluations, which has never been done in GOTM. 60 The remainder of this paper is organized as follows. Development of incorporating CVMix and a Stokes drift module in GOTM is described in Section 2, together with an introduction of a streamlined workflow with Python (python.org) and Jupyter Notebook (jupyter.org) enabled by the newly developed and more flexible configuration capability in GOTM. Evaluation against available theories and observations in four test cases is presented in Section 3. This paper ends with a brief discussion and main conclusions in Section 4. 65 2 Extending the functionality of GOTM

CVMix in GOTM
The Community Vertical Mixing Project (CVMix, Griffies et al., 2015) is a portable vertical mixing software package providing an extensible framework for the development of first-order turbulence closures. It provides a set of subroutines allowing flexible implementation of surface and interior turbulence closures in the K-profile parameterization (KPP, Large et al., 1994;70 Van Roekel et al., 2018) in an ocean general circulation model (OGCM). These subroutines can be assembled in different ways to accommodate different needs (e.g., loop structure, available mean field variables) of the host OGCMs. Recently, a few Langmuir turbulence parameterization schemes based on KPP have also been included in CVMix (e.g., Li et al., 2016;Reichl et al., 2016;. This allows easy implementations and testing of these Langmuir turbulence parameterization schemes in other models. 75 The General Ocean Turbulence Model (GOTM, Umlauf and Burchard, 2005;Umlauf et al., 2014, see updated version on gotm.net) provides a collection of various turbulence closure schemes for the vertical mixing in the ocean and lakes, in particular second-order turbulence closures. The procedures and variables representing the implementation of these models are encapsulated in GOTM's FORTRAN module turbulence as illustrated in Fig. 1. This module can be easily integrated in any existing library structure of a third-party OGCM or lake model. In addition, GOTM can also be used as a stand-alone 80 one-dimensional water column model with flexible configurations to study the hydrodynamic and thermodynamic processes related to vertical mixing in natural waters. In this case, the turbulence routines are called from inside GOTM's main time loop implemented in the central FORTRAN module gotm (see Fig. 1). For these reasons, GOTM provides a useful platform for developing and comparing different ocean and lake vertical mixing parameterizations in both idealized and realistic scenarios (e.g., Burchard and Bolding, 2001;Umlauf and Burchard, 2005;Burchard et al., 2008;Li et al., 2019). 85 As a first step towards extending the capability of GOTM to include Langmuir turbulence parameterizations, CVMix is incorporated in GOTM as an external library, including the above-mentioned first-order closure models of Langmuir turbulence based on KPP (see Fig. 1). Methodically, this is similar to the approach taken in Li et al. (2019) that enables the cross-comparison among a set of Langmuir turbulence parameterization schemes in the single-column setup. Here, a user-interface is developed and described that ensures consistency with the other modules of the GOTM code, such as the mean flow and 90 the meteorological forcing modules (interface module gotm_cvmix, see Fig. 1). Variables passed to CVMix through this interface module include mean flow variables such as the currents, temperature and salinity, surface forcing variables such as the surface friction velocity and surface buoyancy flux, and wave forcing variables such as the Langmuir number and Langmuir enhancement factor. In return, turbulence variables such as the turbulent diffusivity and viscosity are passed back to the GOTM main time loop. 95 Using the new interface to CVMix, it is now possible to directly call the CVMix subroutines, thereby making a range of recent Langmuir turbulence parameterizations directly available in GOTM, and in OGCMs or lake models that use GOTM.
This allows to objectively compare state-of-the-art versions of the KPP model and second moment closure models. Note that the CVMix interface is separate from the GOTM turbulence module (Fig. 1), so additional modifications to the source code of the host model are still needed. However, such modifications are significantly less than would be needed if CVMix were to be 100 directly implemented in the host model.

Stokes drift in GOTM
Stokes drift (see a recent review by van den Bremer and Breivik, 2018, and references therein) is a key property of ocean surface waves that is crucial for the dynamics and parameterizations of Langmuir turbulence (e.g., McWilliams et al., 1997;Li et al., 2019). To provide the necessary information of the Stokes drift for various Langmuir turbulence parameterizations, 105 Stokes drift variables and a few options to configure the Stokes drift are implemented as a new module in GOTM. These different options of Stokes drift provide a way to test the sensitivity of a Langmuir parameterization to the uncertainties in the estimate of Stokes drift.
The most flexible option is to directly read in the Stokes drift profiles from a file (see FORTRAN module stokes_drift in Fig. 1). The Stokes drift profiles can be either computed from the wave spectrum of direct measurements and wave models, 110 or estimated from some empirical relations. Tools for generating the Stokes drift input file for GOTM from various sources are provided on Github (github.com/qingli411/gotmtool). To assist the development and testing of Langmuir turbulence parameterizations, two idealized options are also implemented (Fig. 1). The first option assumes a monochromatic surface wave, for which the Stokes drift is an exponentially decaying profile with depth defined by the surface value u S 0 and a decay depth scale δ S , where z ≤ 0 is the water depth. The exponential profiles have been used in many idealized large eddy simulations of Langmuir turbulence (e.g., McWilliams et al., 1997;Grant and Belcher, 2009). The second option assumes a Stokes drift profile that depends only on the wind, derived from a set of empirical relations and assumptions (the "theory wave" approach of Li et al., 2017, see their Eq. (25)). This "theory wave" approach estimates the Stokes drift profile from the wind assuming a f −5 (where 120 f is the frequency) spectral shape (sea also Breivik et al., 2016) with the directional spread correction of Webb and Fox-Kemper (2015). The surface value and integrated transport of Stokes drift are estimated using empirical relations. By estimating the Stokes drift from the local wind, the contribution of swell is not explicitly represented, except a constant magnitude loss coefficient tuned against a WAVEWATCH III global wave hindcast simulation to represent the reduction effect by swell (Webb and Fox-Kemper, 2015). It has been shown to provide enough information of the ocean surface waves to allow a reasonable 125 representation of the effects of Langmuir turbulence in OGCMs without coupling with a wave model through either a KPP variant  or an energetics based planetary boundary layer scheme (Reichl and Hallberg, 2018;Reichl and Li, 2019). For both the exponential and the "theory wave" options, the controlling parameters (surface value and decay depth of Stokes drift in the former case and surface wind in the latter) can be set to constant or read from a file.
It should be noted that in this study the Stokes drift profile is only used in the Langmuir turbulence parameterizations without 130 being integrated into the mean flow module in GOTM (e.g., Coriolis-Stokes force). This is consistent with the requirement of the two KPP-based Langmuir turbulence parameterizations used here. For other Langmuir turbulence parameterizations, direct modifications of the mean flow equations in GOTM due to Stokes drift may be necessary (see, e.g., Appendix A of Li et al., 2019), which is straightforward as the full Stokes drift is now readily available in GOTM.

135
GOTM now supports the human-readable data-serialization language YAML (yaml.org) for the configuration of parameters.
YAML is both more user-friendly and more developer-friendly than the FORTRAN namelist originally used in GOTM for the configuration. It has a clean and minimal syntax and is easy for extension and maintenance. The default values and the documentation of the configuration parameters of GOTM are now stored in the source code, so that the configurations can be generated from the compiled GOTM executable. This eliminates the need to save the default values of configuration parameters 140 in a separate file, e.g., an XML file (w3.org/TR/xml11), as previously used in GOTM as well as many other ocean models. This is a useful feature especially for the development and maintenance of GOTM, since the configurations are always consistent with the source code. For GOTM users, this also guarantees that a compatible configuration file is always available whenever the source code is updated, which would require extra efforts if FORTRAN namelist files were used.
The Generic Length Scale (GLS, Umlauf and Burchard, 2003) scheme in the k-ε formulation with the weak-equilibrium 160 stability function by Canuto et al. (2001), using a steady-state Richardson number of Ri st = 0.25, denoted as GLS-C01A hereafter, is used as a reference. The parameters used in GLS-C01A are summarized in Table 2.   In the following sections, GOTM simulations with the above four vertical mixing schemes will be compared with available theoretical scalings or observations in four different test cases. We use a Cartesian coordinate system with x and y denoting the horizontal coordinates, z the vertical (upward) coordinate, and u, v and w the corresponding components of the velocity.

Idealized entrainment
The first test case is an idealized wind stress-driven entrainment case with no rotation, in which the OSBL gradually entrains into an underlying non-turbulent region with constant stable stratification. The GOTM simulation results can be directly compared with the relation derived from laboratory experiments (e.g., Price, 1979), in which the time evolution of the mixed layer where R v ≈ 0.6 is the bulk Richardson number, u * the water side surface friction velocity, and N 0 the initial buoyancy frequency. See more discussion on this relation and the model configuration in Umlauf and Burchard (2005). Earth's rotation is not considered.
To test the effect of parameterizing the Langmuir turbulence enhanced mixing and the Langmuir turbulence enhanced entrainment in KPPLT-VR12 and KPPLT-LF17, we assume a Stokes drift in the wind direction (here in the x-direction) that ex-180 ponentially decays with depth following Eq. (1), with a surface value of u S 0 = |u S 0 | = 0.11 m s −1 and decay scale of δ S = 5 m. This corresponds to a turbulent Langmuir number, at which Langmuir turbulence has a prominent influence on the turbulent mixing in the mixed layer (McWilliams et al., 1997).

Ocean Station Papa
The meteorological and oceanic observations at the Ocean Station Papa (OSPapa; 50.1 • N, 144.9 • W; pmel.noaa.gov/OCS/Papa) have been used to evaluate the performance of OSBL turbulent mixing schemes in many studies (e.g., Martin, 1985;Large et al., 1994;Kantha and Clayson, 1994;D'Alessio et al., 1998;Burchard and Bolding, 2001;Acreman and Jeffery, 2007), focusing 210 mostly on the year 1961. Recent measurements of ocean surface waves at OSPapa (Thomson et al., 2013) allow us to evaluate the effects of Langmuir turbulence parameterizations and assess the importance of Langmuir turbulence at this site.
Here we use the temperature and salinity mooring data at OSPapa to initialize the GOTM simulations from rest in a 150 m vertical domain with 150 vertical grid cells. The time step is 60 s. Surface boundary conditions are set by the hourly surface fluxes data from March 21, 2012 to March 20, 2013. Throughout the year the Jerlov water type II (Paulson and Simpson, 1977) 215 is assumed. Half-hourly wave spectral data collected using the Datawell Waverider buoy (cdip.ucsd.edu/metadata/166p1), and binned into n = 64 frequency bands with f 1 = 0.025 Hz and f 64 = 0.58 Hz, is used to estimate the Stokes drift for the Langmuir turbulence parameterizations. The Stokes drift profile is estimated from the band wave energy density spectrum S i according to:

220
where f i is the band center frequency, ∆f i the bandwidth,ê W i a unit vector in the band mean direction and g the gravity acceleration. The grid cell averaged value is computed following Appendix B of Harcourt and D'Asaro (2008). For simplicity, we are ignoring the effect of wave spreading, which may lead to :::: result :: in : an overestimation of the Stokes drift (e.g., Webb and Fox-Kemper, 2015). Note that unlike Li et al. (2019), we are not including the contribution of a f −5 spectral tail beyond the cutoff frequency. A spectral tail contributes more to the surface value of Stokes drift, but much less to the surface layer

230
The annually-averaged net heat flux and freshwater flux over this one-year period are 31.7 W m −2 and 12.9 mg m −2 s −1 , respectively. Such imbalance of the heat and freshwater fluxes would increase the temperature of a water column of 100 m by about 2.5 • C over a year, and decrease the salinity by about 0.1 g kg −1 . 1 To directly compare with the temperature and salinity measurements in single-column simulations over a long period, such imbalance in the heat and freshwater fluxes needs to be compensated by careful adjustments to account for the effects of vertical advection and lateral processes (e.g., Large, 235 1996). Our focus here is to show the effects of parameterizing Langmuir turbulence in GOTM. Therefore, instead of trying to balance the heat and freshwater budget in a rather empirical way, we break the seasonal cycle into four relatively shorter stages (see Fig. 4). These four stages roughly represent (I) the spring restratification, (II) stable forcing in summer, (III) mixed layer entrainment in fall and winter, and (IV) preconditioning for restratification in winter, each of the stages being initialized by observed temperature and salinity profiles. In this way, the differences between different vertical mixing schemes can be 240 shown under different forcing regimes, using the observation as a reference, and the accumulative effects of the ignored vertical advection and lateral processes are reduced.
As shown in Fig. 4e, all four vertical mixing schemes predict warmer SST than the observation throughout the year, especially in stages II and III, and saltier SSS in stage III and slightly fresher SST in stage IV. This is likely a result of missing the vertical advection and lateral processes. Correspondingly, the MLD in all four vertical mixing schemes except KPPLT-LF17 is VR12 and KPPLT-LF17 are strongest during mixed layer deepening in the fall (stage III) and the sporadic mixing events when the mixed layer is shallow (stage II). Its effects are weaker when the mixed layer is deep in winter (stage I and IV), even though both the winds and waves are stronger in winter (Fig. 4c,d). Such effects of Langmuir turbulence in KPP-based 265 parameterizations are significant as compared to the difference between KPP-CVMix and GLS-C01A.
To test the sensitivity of the Langmuir turbulence parameterizations to the uncertainties in estimating the Stokes drift, we repeat the GOTM simulations with Stokes drift estimated from (i) the "theory wave" approach in , and (ii) an idealized exponential profile assuming δ S = 5 m and La t = 0.3 in Eqs. (1) and (3). The "theory wave" estimate of Stokes drift is often an underestimate as the effects of swell are largely ignored. The exponential profile, which represents an idealized 270 swell (but in the direction of the local wind), is likely an overestimation of the Stokes drift for most cases, especially given that the Stokes drift in the real ocean typically decay much faster than exponential Fox-Kemper, 2011, 2015). This is indeed the case for OSPapa, as shown in Fig. 7 which compares the distributions of the Stokes drift profiles estimated from the "theory wave" approach (blue) and an exponential profile (green) with that computed from the observed wave spectrum (see Fig. 4c,d). Note that there is a seasonal variation of the relation between Stokes drift and the wind at OSPapa that is better captured by the "theory wave" estimate but not in the exponential profile estimate (stage II).
Figs. 8a,b compare the time series of La −2 t and La −2 SL from the two estimates of Stokes drift with that computed from the observed wave spectrum, where La SL is the surface layer averaged Langmuir number (Harcourt and D'Asaro, 2008),

280
in which u S SL is the surface layer averaged Stokes drift and u S ref is a reference Stokes drift at the base of the mixed layer. The "theory wave" approach of  gives a slightly smaller estimate of the Stokes drift and the exponential profile approach gives a much bigger estimate for most of the times, especially for the surface layer averaged values. Panels (c)-(e) compare the simulated SST, SSS and MLD between the three simulations. It is seen that KPPLT-VR12 and KPPLT-LF17 give very similar results with both the two estimates of Stokes drift and that computed from the observed wave spectrum, except the 285 SST in stage II where the exponential profile approach significantly overestimates the Stokes drift, and therefore Langmuirenhanced mixing and near-surface cooling. This is consistent with the findings of , suggesting that we may use the "theory wave" estimate of Stokes drift for KPPLT-VR12 and KPPLT-LF17 in GOTM simulations for cases where sufficient wave measurements are not available, such as for the two cases discussed in the following sections. It should be noted that both KPPLT-VR12 and KPPLT-LF17 use only the Langmuir numbers to parameterize the effects of Langmuir turbulence, which are 290 relatively insensitive to the exact profile of Stokes drift. For vertical mixing schemes that depends on the full profile of Stokes drift (e.g., Harcourt, 2013Harcourt, , 2015, this "theory wave" estimate might not be sufficient.

FLEX
The Fladen Ground Experiment (FLEX) test case is based on an intensive field campaign carried out in spring 1976 in the northern North Sea at 58 • 55'N and 0 • 32'E at a depth of about 145 m. Between April 6 and June 13 regular CTD (conductivity-295 temperature-depth) profiles were sampled that were compiled by Soetje and Huber (1980) into vertical profiles of potential temperature and salinity. These profiles show the transition from fully mixed to stratified conditions in the upper half of the water column with a top to mid-depth temperature difference of 4 K. The salinity stratification remains weak with a maximum bottom to top difference of 0.1 g kg −1 . Bottom-generated turbulence due to weak tidal currents keeps the lower half of the water column mixed and supports deepening of the thermocline. The meteorological forcing was highly variable including 300 several storms that led to intermittent mixed layer deepening. Ship-based meteorological data for wind speed, dry and wet air temperature, air pressure, short-wave radiation and long-wave back radiation are available. Since lateral advection is weak in this region and the development of thermal stratification depends on a subtle balance of stratifying forces of surface warming and de-stratifying forces of wind and tidal mixing, the FLEX data set has become a standard test case for surface mixed layer models (e.g., Friedrich, 1983;Frey, 1991;Burchard and Baumert, 1995;Burchard et al., 2006) including GOTM.

305
Since the focus of the present model development is on the surface boundary layer and CVMix does not contain a bottom boundary layer module at the moment, we simulate the FLEX test case with our four vertical mixing schemes without tidal forcing. However, to illustrate the relative effect of Langmuir turbulence versus tidal forcing, we also run GLS-C01A with tidal forcing as a reference.
Here all five GOTM simulations are initialized with temperature and salinity profiles from April 6, 1976, and run through 310 Jun 7, 1976, forced by hourly meteorological data. The surface heat flux is computed internally in GOTM from the meteorological data following Fairall et al. (1996). The surface freshwater flux is ignored. The local depth is 145 m, resolved with 145 evenly distributed grid cells. A time step of 360 s is used, and 3-hourly output is analyzed here. Since stratification is dominated by temperature in this case, the salinity field is relaxed towards the observations on a time scale of 48 hours, and we focus our discussion on the temperature field.

315
There are no direct measurements of ocean surface waves in the FLEX case, therefore we are using the "theory wave" approximation of  to estimate the Stokes drift from the wind speed, assuming that wind and waves are aligned.
As demonstrated in the previous section, this approximation provides a reasonable estimate of the Stokes drift in our Langmuir turbulence parameterizations KPPLT-VR12 and KPPLT-LF17.

Gotland Basin
Finally, we tested the four vertical turbulent mixing schemes in the central eastern Gotland Basin of the Baltic Sea (57.3 • N, 20.0 • E). The primary goal of this analysis is to test the performance of these schemes in a multi-year simulation of a non-tidal 350 basin that remains stratified throughout the year due to a permanent halocline (Feistel et al., 2008). A detailed description of the GOTM setup for the Gotland Basin is described in Burchard et al. (2006).

390
For periods outside the winter months, however, all four vertical mixing schemes give quite similar results, especially during spring (Fig. 12b,f) and summer (Fig. 12c,g). The difference from the observed temperature profiles may be related to processes missing in the single column model used in our study, for example those associated with the basin-scale doming of the density structure (see, e.g., Fig. 9  OSPapa case (Fig. 6), and the increased MLDs discussed in the context of Fig. 12d,h above, the most apparent differences among the four schemes are found in the phase of fall-to-winter mixed layer deepening, with KPPLT-LF17 giving the strongest deepening and KPP-CVMix the weakest. This is expected from both the design of the schemes and from results of previous tests. Overall the three KPP variants in the newly implemented CVMix module perform reasonably well in this relative long simulation.  (not shown). Note, however, that the relative effects of including Langmuir turbulence in KPPLT-VR12 and KPPLT-LF17 as shown in Fig. 6 do not change with the vertical resolution.
In comparison, as shown in Fig. 15, GLS-C01A is much less sensitive to the vertical resolution than KPP-CVMix. Fig. 15d   415 shows that during the phases of strong mixed layer deepening (in particular, phase III), differences in MLD and mixed-layer temperatures are approximately an order of magnitude smaller compared to KPP-CVMix if the grid spacing is increased from 1 m to 5 m. However, GLS-C01A is noticeably more sensitive to the time step than KPP-CVMix (Fig. 15b,c). This is likely related to the fact that GLS-C01A solves prognostic equations for the turbulent velocity scale and length scale, as compared to the diagnostic algorithms in KPP. However, the numerical solution of partial differential equations in GLS-C01A has the 420 important advantage that numerical errors are guaranteed to decrease according to the well-defined convergence properties of the numerical schemes if the vertical grid spacing is reduced. For the algorithms in KPP-CVMix, this property cannot be proved.
We also compared the execution (CPU) times in GLS-C01A and KPP-CVMix for all the relevant subroutines in each scheme, excluding common subroutines such as the mean flow equations and data input/output. Surprisingly, the CPU time required for 425 a single time step in KPP-CVMix turned out to be 3-4 times larger than for the second-order turbulence model in GLS-C01A. It should be noted that this relative timing information only applies for the specific numerical implementations in GOTM and CVMix used in our study. The actual performance of the corresponding second-order and KPP mixing schemes may be improved by optimizing the loop structure and the time step, which may also change the relative performance of these models.
From a practical point of view, the execution times will also depend on the time step and grid spacing chosen to yield a desired 430 accuracy. As shown above, the second-moment closures in GOTM will tolerate the use of coarser vertical grids, whereas KPP-CVMix will provide sufficient accuracy also for larger time steps. Overall, therefore, our study does not show any clear advantage in computational costs for either model. This is a rather remarkable result as KPP-type models are generally believed to be more robust, and therefore preferable, in coarse-resolution global modeling. based configuration lists, as well as for a straightforward integration of additional Langmuir turbulence parameterizations in the future. For example, future development of CVMix, which has been widely used in many ocean climate models, can now be easily incorporated in GOTM through the new CVMix module, and are therefore available to other ocean models that use GOTM as their turbulence library with modest code changes. In addition, even though we have demonstrated in Section 3.2 that the two variants of KPP with Langmuir turbulence are not particularly sensitive to the details of the Stokes drift profiles 445 (consistent with previous studies), the full Stokes drift profiles are provided by the Stokes drift module, which will facilitate future development and incorporation of Langmuir turbulence parameterizations in the GOTM framework, such as the second moment closures of Langmuir turbulence by Harcourt (2013Harcourt ( , 2015. Using these two new modules in GOTM, three variants of KPP in CVMix, with and without Langmuir turbulence, were tested and compared with a second-moment turbulence closure scheme based on the GLS-framework in the k-ε formulation. We investigated four test cases, and evaluated the model performance against available theoretical scalings and observations. These four test cases included an idealized wind-driven entrainment case and three more realistic cases: the Ocean Station Papa, the northern North Sea, and the central Gotland Sea, each focusing on slightly different aspects of the vertical mixing processes. The results are consistent with previous studies of Langmuir turbulence effects in KPP (e.g., Li et al., 2019). Although the degree to which we can evaluate these vertical mixing schemes is still limited by the use 455 of single column simulations, interesting conclusions can still be drawn from the direct comparison with available theoretical scalings and observations, which has never been done particularly for KPPLT-VR12 and KPPLT-LF17. The effects of Langmuir turbulence as represented in KPPLT-VR12 and KPPLT-LF17 in these test cases are most important during periods when the mixed layer is deepening. Such effects on the simulated SST, SSS and MLD in a single column setup can sometimes be comparable to the effects of the missing advection and lateral processes at Ocean Station Papa (inferred from the mismatch 460 between modeled and observed). The magnitude of such effects also appears to be comparable to or even larger than the effects of tidal forcing in the northern North Sea in the FLEX case. These results highlight the importance of correctly representing the effects of Langmuir turbulence in an ocean vertical mixing scheme. Compared with the k-ε scheme, all three variants of KPP suffer from much higher sensitivity to the vertical resolution, consistent with previous studies (e.g., Van Roekel et al., 2018;Li et al., 2019), but lower sensitivity to the time step.

465
The newly developed YAML-based configuration in GOTM also enables an easier GOTM workflow using Python and Jupyter Notebook, both gaining popularity rapidly over the years in the broader scientific community. This provides an interface to use GOTM interactively in the Jupyter Notebook environment, especially as a single water column model which has been extremely useful in both parameterization development and evaluation (e.g., Burchard and Bolding, 2001;Burchard et al., 2008;Reichl et al., 2016;Li et al., 2019) and parameter space exploration (e.g., Li et al., 2019;Dong et al., 2020). As an 470 example, all the simulations and figures in this paper are conducted and produced in the Jupyter Notebook environment using this new workflow. Both the source code for this new workflow and the Jupyter notebooks for this paper are available online for maximum reproducibility (see the Code and data availability below for more details).
This study represents the initial steps extending the capability of GOTM to include Langmuir turbulence parameterizations, although we note that GOTM has already been used in previous studies in developing and evaluating Langmuir turbulence 475 parameterizations (e.g., Reichl et al., 2016;Li et al., 2019). Even though only a limited set of KPP variants is included, the development here facilitates future incorporation of second moment closure schemes of Langmuir turbulence (e.g., Harcourt, 2013Harcourt, , 2015, as well as the development and evaluation of new Langmuir turbulence parameterizations in the GOTM framework. The Python and Jupyter Notebook based GOTM workflow enabled by the YAML-based configuration also makes future applications of GOTM easier.

480
It should be noted that GOTM can also be coupled to the frequently used Framework for Aquatic Biogeochemical Models (FABM, Bruggeman and Bolding, 2014) to study the evolution of marine ecosystems. The effect of different turbulence models for ecosystem-related questions can now be evaluated within a single modeling framework. However, a systematic evaluation of such effect is beyond the scope of this work and is left for future study. ("Energy-consistent atmosphere ocean coupling"), embedded in the Collaborative Research Centre TRR 181 "Energy Transfers in Atmosphere and Ocean". KK and HB are supported by TRR 181 as well. This project did also fund parts of the integration of CVMix into GOTM carried out by KB and JB who recieved further support from IOW to improve the YAML-based GOTM configurations.