Reanalysis of vertical mixing in mesocosm experiments: PeECE III and KOSMOS 2013

Controlled manipulation of environmental conditions within large enclosures in the ocean, so-called pelagic mesocosms, has become a standard method to explore potential responses of marine plankton communities to anthropogenic change. Among the challenges of interpreting mesocosm data is the often uncertain role of vertical mixing, which usually is not observed directly. To account for mixing nonetheless, two pragmatic assumptions are common: either that the water column is homogeneously mixed or that it is divided into two water bodies with a horizontal barrier inhibiting turbulent exchange. In this study, we present a model-based reanalysis of vertical turbulent diffusion in the mesocosm experiments PeECE III and KOSMOS 2013. Our diffusivity estimates indicate intermittent mixing events along with stagnating periods and yield simulated temperature and salinity profiles that are consistent with the observations. Here, we provide the respective diffusivities as a comprehensive data product in the Network Common Data Format (NetCDF). This data product will help to guide forthcoming model studies that aim at deepening our understanding of biogeochemical processes in the PeECE III and KOSMOS 2013 mesocosms, such as the CO₂-related changes in marine carbon export. In addition, we make our model code available, providing an adjustable tool to simulate vertical mixing in any other pelagic mesocosm. The data product and the model code are available at https://doi.org/10.1594/PANGAEA.905311 (Mathesius et al., 2019).

1 Introduction

In a rapidly changing world (Steffen et al., 2011, 2018; Stocker et al., 2014), a better understanding of anthropogenic pressures on marine ecosystems and biogeochemical cycles is needed, so that major risks can be identified in advance and mitigated. The response of marine plankton communities to anthropogenic change is often explored in so-called pelagic mesocosms, large polyethylene bags, either drifting or mounted in coastal waters, usually 2 m in diameter and up to 20 m long. They are filled with local seawater and contain a natural plankton community (Riebesell et al., 2008, 2013). Mesocosms are big enough to prevent many “bottle effects” that are existent in laboratory experiments, while at the same time they are small enough to allow for a cost-efficient controlled manipulation of environmental conditions and frequent monitoring of species composition and biogeochemistry. Therefore, they provide unique insights into possible responses of planktonic ecosystems to anthropogenic pressures that cannot be obtained from laboratory studies. Past mesocosm experiments investigated, for example, the effects of decreasing pH (Riebesell et al., 2008, 2013; Gazeau et al., 2016; Engel et al., 2014; Galgani et al., 2014; Archer et al., 2018), nutrient increase (Schwier et al., 2017; Micheli, 1999), and warming (Lewandowska et al., 2014; Sommer et al., 2015; Wohlers et al., 2009). Ocean acidification has been found to have a significant impact on organic matter stoichiometry (Riebesell et al., 2007; Schulz et al., 2008), plankton community composition (Bach et al., 2016; Riebesell et al., 2017; Boxhammer et al., 2018), primary production, and carbon export (Egge et al., 2009; Riebesell et al., 2007). Even though mesocosms provide a well-controlled environment, where causes and consequences are much easier to disentangle than in the field, uncertainties remain. The observed changes in biogeochemistry and community structure are usually a result of numerous unobserved processes, so there is often room for interpretation, with more than one plausible hypothesis regarding the key events that took place during the experiment. Models are a valuable complement to statistical data analyses, because they can explicitly resolve and investigate unobserved processes and therefore test different assumptions and hypotheses. In this study, we focus on the uncertain role of vertical mixing, which is very hard or even impossible to measure directly without effectively mixing the water column by introducing respective measurement devices. So far, the interpretation of some of the most influential mesocosm experiments could not take vertical mixing into account, since data on mixing were not available. However, vertical mixing can be an important factor that affects sinking and resuspension of particulate matter, dilution and entrainment of nutrients, and air–sea CO₂ fluxes. By distributing phytoplankton to different depths, vertical mixing also influences the light-exposure and consequently light-sensitive physiological processes of the phytoplankton, such as photo-acclimation, which has an effect on primary production. Accounting for vertical mixing in the interpretation of mesocosm data might provide a better understanding of the observed biogeochemical changes. To this end, we developed a one-dimensional mesocosm model that is able to retrace or reanalyze vertical turbulent diffusivities based on temperature, salinity, and solar radiation data. These data sets are available for nearly all mesocosm experiments, so that this model can be used to simulate vertical mixing in almost any pelagic mesocosm. In this paper, we present simulations of the mesocosm experiments PeECE III (Pelagic Ecosystem CO₂ Enrichment III; Schulz et al., 2008) and KOSMOS 2013 (Kiel Off-Shore Mesocosms for Future Ocean Simulations 2013; Bach et al., 2016). These experiments exposed marine plankton communities to elevated CO₂ levels and provided highly valuable insights into the plankton community response to ocean acidification. There are indications that the vertical mixing patterns in the mesocosm experiments PeECE III and KOSMOS 2013 are complex, since a variety of forcing mechanisms are at work. Among them are (1) thermal stratification due to solar radiation; (2) a diurnal cycle in air–sea heat and associated buoyancy fluxes, where surface waters are cooled at night, inducing destabilization of the water column and convection; and (3) wind friction driving shear and waves which impinge on the flexible mesocosm walls, thereby turbulently pushing water up and down. Despite evidence that substantial mixing as well as strong stratification have occurred in the mesocosms, quantitative estimates have not been available. Here, we publish the model code and our estimated turbulent diffusivities for the PeECE III and KOSMOS 2013 mesocosm experiments, available at https://doi.org/10.1594/PANGAEA.905311 (Mathesius et al., 2019). Furthermore, we demonstrate that these estimates provide realistic profiles of temperature and salinity that are consistent with the observations. The mixing model can easily be coupled to any plankton ecosystem mesocosm model (e.g., the model of Krishna and Schartau, 2017). In a follow-up study, we will present the effect of vertical mixing on biogeochemical tracers in mesocosms.

2 Methods

2.1 Data

In this paper, we show the simulations of two mesocosm experiments, PeECE III (Schulz et al., 2008; Riebesell et al., 2007) and KOSMOS 2013 (Bach et al., 2016). The experiments differ in duration, length of mesocosms, environmental conditions, and initial conditions. At the beginning, the PeECE III mesocosms are strongly stratified, while the KOSMOS 2013 mesocosms are fully mixed and then become partially stratified during the course of the experiment. Thus, these two mesocosm experiments cover a spectrum of different mixing patterns and provide a good test ground for our mixing model.

The PeECE III experiment was conducted in the fjord of Bergen, Norway, from 5 May to 15 June 2005. Nine mesocosms of 8.5 m length were deployed in the fjord and sampled every day. Temperature and salinity profiles were measured daily by a conductivity–temperature–depth (CTD) sensor. Before the experiment started, fresh water was mixed into the upper 5.5 m to create an artificial stratification. Throughout the experiment the upper surface layer was mixed by a pump to guarantee a homogeneous distribution of dissolved compounds within the mixed layer. Solar radiation data were provided by the Geophysical Institute of the University of Bergen (Olseth et al., 2006). Data on chlorophyll a concentrations were provided via PANGAEA by the PeECE III team (2008).

For the KOSMOS 2013 experiment, 10 mesocosms of 19 m length were deployed in the Gullmar Fjord, located near Kristineberg at the west coast of Sweden, from 7 March to 28 June 2013. CTD measurements and depth-integrated water samples were taken every other day (data provided by Boxhammer et al., 2017). There was no artificial stratification; instead the water column was fully mixed at the start of the experiment. For our simulations, we used solar radiation data from the nearby location of Kristineberg, provided by the Sven Lovén Centre for Marine Sciences of the University of Gothenburg (http://www.weather.loven.gu.se/kristineberg/en/data.shtml, last access: 9 February 2016). Data on chlorophyll a concentrations were provided via PANGAEA by Boxhammer et al. (2017).

For our simulations and model–data comparisons, we linearly interpolate the CTD data to hourly values. Data gaps are filled by linear interpolation as well. Since temperatures were measured only at daytime, the diurnal cycle of temperature was not observed and nightly cooling cannot be resolved in our model (implications are discussed in Sect. 4.5).

2.2 Model description

The model has a vertical resolution of 0.5 m. The temporal resolution is process-specific: 1 s for the turbulent diffusion, 10 min for the convective adjustment, and 1 h for the radiative heating and the surface corrections (see below). As initial profiles for temperature and salinity in each mesocosm, we take the average of the measured temperature and salinity profiles of the first 3 experiment days.

To simulate the temporal evolution of temperature and salinity within the mesocosms, we set up a one-dimensional turbulent diffusion model for these two tracers. We start with the general diffusion equation for a tracer ϕ, which is given by

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(k_z\cdot {\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right),\end{array}$$

where k_z denotes the depth-dependent eddy diffusivity and z the depth. The diffusion equation is solved based on a code of Burkardt (2009) that uses second-order central differences to approximate the second derivative in space and an implicit Euler approximation for the first derivative in time.

Vertical mixing in mesocosms can be influenced by different degrees of stratification of the water column, where the time and depth of stratification (or lack thereof) can be highly variable. Based on the density-dependent diffusivity parametrization of Osborn (1980) that inherently takes stratification-induced inhibition of vertical mixing into account, we parameterize the diffusivity k_z as

$$\begin{array}{}\text{(2)}& k_z={\displaystyle \frac{c}{N_z^{\mathrm{2}}}},\end{array}$$

where the parameter c corresponds to the product of dissipation rate and mixing efficiency of Osborn's original parametrization, and $N_z^{\mathrm{2}}$ is the buoyancy frequency given by

$$\begin{array}{}\text{(3)}& N_z^{\mathrm{2}}={\displaystyle \frac{g}{{\mathit{\rho }}_z}}\cdot {\displaystyle \frac{\partial {\mathit{\rho }}_z}{\partial z}},\end{array}$$

where g denotes the gravitational acceleration and ρ_z the potential density of depth z (with z defined positive downward). Both k_z and $N_z^{\mathrm{2}}$ are calculated for every time step and depth level. In general, the dissipation rate as well as the mixing efficiency can significantly vary in time and space. However, our model assumes that c is constant for all depths and the whole time period of the experiment. The constant c, which is optimized for every single mesocosm independently (see Sect. 2.3), can be interpreted as an averaged value within one mesocosm for the entire experiment time. With this simplification, any temporal and vertical variations in the vertical diffusivities k_z are induced by the buoyancy frequency and are inversely proportional to the vertical density gradient. The density gradients and the diffusivities are calculated each hour from the interpolated observed temperature and salinity profiles. By using the observed values instead of simulated temperature and salinity, we avoid follow-up errors that could grow rapidly. Simulated temperature and salinity profiles only become relevant for the diffusivity calculation when the associated simulated density gradient is negative. In this case, we prescribe that convection takes place and we set k_z to a high value of 0.1 m² s⁻¹, which is in general the maximum value that diffusivities can take in our model. We check every 10 min if conditions for convective mixing are present.

Once per hour, we carry out a surface correction of temperature and salinity (Sect. 2.2.2, 2.2.3) and account for heating by solar radiation (Sect. 2.2.1).

2.2.1 Radiative heat flux

Solar radiation penetrates the water column and warms subsurface waters. To calculate the warming related to the radiative heat flux, we use solar radiation observations (see Sect. 2.1) as incoming solar radiation.

Within each box of the model, the radiative heating rate (RHR) is calculated, following Ohlmann (2003), as

$$\begin{array}{}\text{(4)}& {\mathrm{RHR}}_z={\displaystyle \frac{{\mathrm{SI}}_{\mathrm{in}}-{\mathrm{SI}}_{\mathrm{out}}}{c_p\cdot {\mathit{\rho }}_z\cdot \mathrm{\Delta }z}},\end{array}$$

where SI_in denotes the solar radiation entering the box, SI_out the solar radiation leaving the box, Δz the box thickness (0.5 m), and c_p the specific heat capacity of seawater (c_p = 3990 W s kg⁻¹ K ⁻¹). The amount of solar radiation leaving the box is given by

$$\begin{array}{}\text{(5)}& {\mathrm{SI}}_{\mathrm{out}}={\mathrm{SI}}_{\mathrm{in}}\cdot e^{-\mathit{\alpha }\mathrm{\Delta }z},\end{array}$$

and depends on SI_in and the absorption coefficient α, given as

$$\begin{array}{}\text{(6)}& \mathit{\alpha }={\mathit{\alpha }}_{\mathrm{w}}+{\mathit{\alpha }}_{\mathrm{chl}}\cdot \mathrm{chl},\end{array}$$

where α_w is the absorption coefficient of water (α_w=0.04 m⁻¹), α_chl the absorption coefficient of chlorophyll a (α_chl=0.01 m² (mg chl)⁻¹), and chl the chlorophyll a concentration. Here we use the chlorophyll a concentrations measured during the mesocosm experiments. If these were not available, modeled chlorophyll a concentrations could be used.

The product of RHR_z and Δt is the radiative heating per simulation time step and is added to the respective temperature at depth z.

2.2.2 Temperature surface correction

Apart from the radiative heat flux, the water temperature is strongly influenced by sensible and latent heat fluxes at the sea surface. In our model setup we assume that all changes in total heat content within a mesocosm are caused by the combination of radiative heat fluxes and sea-surface heat fluxes.

To account for the sea-surface heat fluxes, we determine the hourly changes in simulated heat content (Eq. 7) and hourly interpolated observed heat content. The difference in heat content change is converted back to temperature (Eq. 8) and added to the top box of the simulated mesocosm, thus adding missing heat or subtracting excess heat.

The water column heat content, H in units of joules per square meter, is calculated each hour for observed and simulated heat, according to

$$\begin{array}{}\text{(7)}& H=c_p\cdot \underset{z=\mathrm{0}}{\overset{z_{max}}{\int }}{\mathit{\rho }}_z\cdot T_z\mathrm{d}z,\end{array}$$

where ρ_z and T_z are the potential density and the temperature at depth z, respectively. We call the hourly change in observed heat ΔH_obs and the hourly change in simulated heat ΔH_sim. The surface temperature correction term T_corr, which accounts for both sensible and latent heat fluxes, is then given by

$$\begin{array}{}\text{(8)}& T_{\mathrm{corr}}=\left\{\begin{array}{ll}\frac{(\mathrm{\Delta }{H}_{\mathrm{obs}}-\mathrm{\Delta }{H}_{\mathrm{sim}})}{{\mathit{\rho }}_{\mathrm{1}}\cdot c_p\cdot \mathrm{\Delta }z}& |z\le \mathrm{0.5}\mathrm{m}\\ \mathrm{0}& |z>\mathrm{0.5}\mathrm{m},\end{array}\right.\end{array}$$

where ρ₁ is the potential density of the top box of the model.

The complete temperature equation that includes the turbulent thermal flux, radiative heating, and sea-surface heat flux correction is then given by

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(k_z\cdot {\displaystyle \frac{\partial T}{\partial z}}\right)+{\mathrm{RHR}}_z+{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{T}}_{\mathrm{corr}}}{\partial t}},\end{array}$$

where ${\stackrel{\mathrm{̃}}{T}}_{\mathrm{corr}}$ represents the continuous form of T_corr.

2.2.3 Salinity surface correction

Similar to the heat correction above, we also introduce a correction term for salinity that corrects for sea-surface fluxes that are not explicitly resolved in the model. We assume that all changes in salinity are caused by freshwater fluxes at the surface, driven by evaporation or precipitation. Every hour, we add the hourly change of observed depth-integrated salinity, S_corr, to the top box of the model.

$$\begin{array}{}\text{(10)}& S_{\mathrm{corr}}=\left\{\begin{array}{ll}{\int }_{z=\mathrm{0}}^{z_{max}}\mathrm{\Delta }{S}_{\mathrm{obs}}\mathrm{d}z& |z\le \mathrm{0.5}\mathrm{m}\\ \mathrm{0}& |z>\mathrm{0.5}\mathrm{m},\end{array}\right.\end{array}$$

where ΔS_obs denotes the hourly change in observed salinity.

The complete salinity equation that includes the turbulent salinity flux and the sea-surface freshwater flux correction can be written as

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial S}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left(k_z\cdot {\displaystyle \frac{\partial S}{\partial z}}\right)+{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{S}}_{\mathrm{corr}}}{\partial t}},\end{array}$$

where ${\stackrel{\mathrm{̃}}{S}}_{\mathrm{corr}}$ represents the continuous form of S_corr.

2.3 Optimization of the diffusivity parameter

To find the best diffusivity estimates for each mesocosm, we optimize the parameter c of Eq. (2), by minimizing three different cost functions that either regard (i) temperature alone, (ii) only salinity, or (iii) both tracers together. For the optimization we take only the daily measurements into account, not the interpolated hourly values. Since the measurements were taken approximately at noon, we compare them to a 3 h average of the simulated tracer values around noon.

When optimizing c depending on temperature, the averaged cost function J_T is given by

$$\begin{array}{}\text{(12)}& J_T={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_T^{\mathrm{2}}n}}\sum _{i=\mathrm{1}}^n{\left(T_{{\mathrm{obs}}_i}-T_{{\mathrm{sim}}_i}\right)}^{\mathrm{2}},\end{array}$$

where T_obs and T_sim denote the observed and simulated temperature, respectively; σ_T the standard deviation of observed temperature in the mesocosm; and n the number of data points (i.e., the product of the number of experiment days and depth levels).

When optimizing c depending on salinity, the averaged cost function J_S is given by

$$\begin{array}{}\text{(13)}& J_S={\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\sigma }}_S^{\mathrm{2}}n}}\sum _{i=\mathrm{1}}^n{\left(S_{{\mathrm{obs}}_i}-S_{{\mathrm{sim}}_i}\right)}^{\mathrm{2}},\end{array}$$

where S_obs and S_sim denote the observed and simulated salinity, respectively, and σ_S the standard deviation of observed salinity in the mesocosm.

Furthermore, we optimize the diffusivity parameter c for both temperature and salinity simultaneously. In this case the so-called combined cost function J_TS is given by

$$\begin{array}{}\text{(14)}& J_{\mathrm{TS}}=J_T+J_S.\end{array}$$

The cost functions J_T, J_S, and J_TS are computed for different values of c, ranging from 10^−8.5 to 10^−6.5 in logarithmic equidistant steps of 0.1. The optimal value of c is where the cost function reaches a minimum. In the following, c_T denotes the optimal value for c when the cost function is optimized with respect to temperature only (Eq. 12), c_S when the cost function is optimized with respect to salinity only (Eq. 13), and c_TS when the cost function is optimized with respect to temperature and salinity simultaneously (Eq. 14). The respective costs are given by J_T(c_T), J_S(c_S), and J_TS(c_TS).

3 Results

3.1 PeECE III

First, we optimize the diffusivity parameter c by finding the corresponding minimum cost. Figure 1 shows the cost values resulting from the temperature optimization, the salinity optimization, and the optimization that depends on both temperature and salinity. We find a distinct minimum for each optimization and mesocosm. This minimum determines the optimal value for c that we use to calculate the turbulent eddy diffusivities (Eq. 2). In all mesocosms, the best estimates of c resulting from the temperature optimization (c_T) are significantly larger than the best estimates of the salinity optimization and the optimization that depends on temperature and salinity simultaneously (c_S and c_TS, respectively).

Figure 1 Costs of three cost functions depending on parameter c (Eq. 2), for the experiment PeECE III. Black lines illustrate the cost function J_TS, taking into account the model–data misfit of temperature and salinity. Red lines illustrate the cost function J_T, taking into account only temperature errors; blue lines illustrate the cost function J_S, taking into account only salinity errors. The horizontal bars on the x axis mark the range of all mesocosms' best estimates of c for the three optimizations.

Figure 2 Temporal evolution of temperature, salinity, and diffusivity profiles in one exemplary mesocosm of PeECE III. The first row shows observed temperature and salinity; the following rows show simulated temperature, salinity, and diffusivity for (a) the salinity optimization (second row), (b) the temperature optimization (third row), and (c) the optimization that depends on both temperature and salinity (fourth row). Shown are results for mesocosm no. 8; figures for the other PeECE III mesocosms are provided in the Supplement.

Table 1 Best estimates of parameter c (c_T, c_S, c_TS) and associated combined cost that takes temperature and salinity errors into account, J_TS(c_T), J_TS(c_S), and J_TS(c_TS), for all mesocosms of the experiment PeECE III.

Table 2 Mean and standard deviation of observed temperature (T) and salinity (S) as well as root-mean-square errors (RMSEs) of simulated temperature and salinity, based on the simulation with c_TS as the optimal value for parameter c. The listed values are for all mesocosms of the experiment PeECE III.

Figure 2 shows the observed temperature and salinity of one exemplary mesocosm (no. 8) as well as the simulated temperature, salinity, and diffusivity that result from simulations with the three optimizations' best estimates for c (c_S, c_T, and c_TS).

When optimizing c depending on salinity only (Fig. 2, second row), the resulting salinity profiles are very similar to the observations. In comparison, the deviation of simulated temperatures from observed temperatures is much larger, especially from day 16 onwards. When optimizing c depending on temperature only (Fig. 2, third row), the simulated temperature profiles are very similar to the observations. However, there is some excess heat in the lower 2 m, especially after the surface heat peak. The simulated salinity pattern differs substantially from the observed pattern, starting already during the first 4 d of the simulation with a rapid weakening of the halocline. When temperature and salinity are both taken into account in the optimization (Fig. 2, bottom row), the results for temperature, salinity, and diffusivity are very similar to the results of the salinity optimization. The simulated temperatures deviate more from the observations than in the temperature optimization, whereas the simulated salinities are almost as accurate as in the salinity optimization.

The reason for the higher impact of salinity is that, even though diffusivities resulting from a c value higher than c_TS improve temperature results, they also cause a quick dissolution of the halocline. This introduces large model errors and thus high costs from the salinity term in Eq. (14). In the opposite case, diffusivities resulting from a c value lower than c_TS preserve the halocline well and at the same time lead to a temperature pattern that still resembles the basic characteristics of the observed pattern, so costs from the temperature term in Eq. (14) are relatively small. Thus, c_TS is much closer to c_S than to c_T.

To underpin these findings, we calculate the cost J_TS depending on temperature and salinity for each optimized simulation. Table 1 shows the three optimizations' best estimates for parameter c and the corresponding cost from the combined cost function, J_TS(c_T), J_TS(c_S), and J_TS(c_TS), for each mesocosm. In general, c_S and c_TS are very close to each other, in one mesocosm even identical, whereas all mesocosms' estimated c_T values are substantially higher. The cost values J_TS(c_T) are in all mesocosms larger than the cost values J_TS(c_S). The combined cost resulting from c_TS, J_TS(c_TS), provides the lowest cost. Thus, the model yields the best results when both temperature and salinity are included in the optimization of diffusivities.

The optimized model simulations are further evaluated by a comparison of the root-mean-square error (RMSE). Here we focus on the best fitted simulation resulting from the simultaneous temperature and salinity optimization. For each mesocosm, Table 2 shows mean and standard deviation of observed temperature and salinity as well as the RMSE of simulated temperature and salinity. The mesocosms' mean temperatures range between 9.93 and 9.98 ^∘C with a standard deviation of 0.76–0.82 ^∘C, while the mean salinities range between 31.03 and 31.17 PSU with a standard deviation of 0.41–0.59 PSU. The RMSE of temperature ranges from 0.24 to 0.30 ^∘C, whereas the RMSE of salinity ranges from 0.09 to 0.15 PSU. The highest RMSE of temperature and salinity in all mesocosms of the PeECE III experiment constitutes only 36.6 % and 28.5 % of the standard deviation of temperature and salinity observations, respectively.

3.2 KOSMOS 2013

For the KOSMOS 2013 experiment, we determine the optimal diffusivities by performing three optimizations, depending on either temperature or salinity only or depending on both (see Fig. 3). The resulting cost values for each parameter variation are quite different from the costs of the previous experiment.

In the salinity optimization, the cost minima are not as pronounced as they are in PeECE III. Best estimates of c_S are higher for all KOSMOS 2013 mesocosms, if compared with those of the PeECE III experiment. Furthermore, the range of c_S values is substantially larger. The opposite is the case in the temperature optimization, where the cost minima are very pronounced, the range of c_T is narrower, and the cost curves of all mesocosms are almost identical. The c_TS values are much higher than in PeECE III, thereby partly overlapping with the estimated c_T and c_S values.

Table 3 lists the best estimates for parameter c (c_T, c_S, and c_TS) and the corresponding cost J_TS for all mesocosms. In four mesocosms, c_T equals c_TS; i.e., the temperature optimization provides the same optimal c value and therefore the same diffusivities as the simultaneous temperature and salinity optimization. In two mesocosms, c_S equals c_TS. In the remaining mesocosms, c_TS is between c_T and c_S and provides the lowest cost. When c_TS equals c_S or c_T, the corresponding combined cost values are identical. However, whenever the optimal c values are not identical, the simultaneous temperature and salinity optimization leads to the lowest model–data misfit.

Figure 3 Costs of three cost functions depending on parameter c (Eq. 2), for the experiment KOSMOS 2013. Black lines illustrate the cost function J_TS, taking into account the model–data misfit of temperature and salinity. Red lines illustrate the cost function J_T, taking into account only temperature errors; blue lines illustrate the cost function J_S, taking into account only salinity errors. The horizontal bars on the x axis mark the range of all mesocosms' best estimates of c for the three optimizations.

Figure 4 shows the observed temperature and salinity of one exemplary mesocosm (no. 10) as well as the simulated temperature, salinity, and diffusivity that result from simulations with the three optimizations' best estimates for c. The depicted mesocosm (no. 10) is the one that shows the largest spread between c_S, c_T, and c_TS.

Figure 4 Temporal evolution of temperature, salinity, and diffusivity profiles in one exemplary mesocosm of KOSMOS 2013. The first row shows observed temperature and salinity; the following rows show simulated temperature, salinity, and diffusivity for (a) the salinity optimization (second row), (b) the temperature optimization (third row), and (c) the optimization that depends on both temperature and salinity (fourth row). Shown are results for mesocosm no. 10; figures for the other KOSMOS 2013 mesocosms are provided in the Supplement.

The observed temperature and salinity profiles are substantially different from the conditions found in the PeECE III experiment. Within the first 5 weeks, the whole water column is well mixed and characterized by temperatures below 4.8 ^∘C and salinities below 29.5 PSU. Over the course of the experiment, temperature and salinity increase to up to 16.8 ^∘C and 29.66 PSU, respectively, and a pronounced thermal stratification is established for several weeks.

In all three optimizations, the model is able to reproduce the observed features very well. Notably, the simulated temperature and salinity profiles of the simultaneous temperature and salinity optimization are more similar to the profiles of the temperature optimization, whereas in PeECE III they were closer to the profiles of the salinity optimization.

The diffusivities resulting from the simulations with c_T and c_TS are very similar as well, while the diffusivities of the salinity optimization are significantly larger. The right column of Fig. 4 shows that the simulated water column is well mixed during the first 5 weeks, followed by an increasingly distinct stratification with lower diffusivities between surface and bottom water masses. However, there are periods of deep convective mixing, reaching almost down to the mesocosm bottom. Like in the optimized simulations of PeECE III, heat at the bottom is slightly overestimated in all simulations.

The mesocosms' mean temperatures range between 7.36 and 7.40 ^∘C with a standard deviation of 4.49–4.53 ^∘C, while the mean salinities range between 29.10 and 29.45 PSU with a standard deviation of 0.07–0.11 PSU (Table 4). We further evaluate the optimized model simulations by a comparison of the root-mean-square error (RMSE), focusing on the best-fitted simulations that result from the simultaneous temperature and salinity optimization. The RMSE of temperature ranges between 1.00 and 1.09 ^∘C, which is much higher than in PeECE III. However, this does not mean that the model performance is worse, since the RMSE should be interpreted in the context of the overall variation in observed temperature. Putting the RMSE of temperature into perspective, it constitutes only 24.2 % of the variability (standard deviation) of the observed temperature, which is lower than the corresponding ratio in the PeECE III mesocosms (see Sect. 3.1). In contrast, the RMSE of salinity is lower than in PeECE III, ranging from 0.02 to 0.04 PSU, while the ratio of the highest RMSE and the standard deviation of the observed salinity is 40.6 % and thus larger than in PeECE III.

Table 3 Best estimates of parameter c and associated combined cost that takes temperature and salinity errors into account, J_TS(c_T), J_TS(c_S), and J_TS(c_TS), for all mesocosms of the experiment KOSMOS 2013.

Table 4 Mean and standard deviation of observed temperature (T) and salinity (S) as well as root-mean-square errors (RMSEs) of simulated temperature and salinity, based on the simulation with c_TS as the optimal value for parameter c. The listed values are for all mesocosms of the experiment KOSMOS 2013.

4 Discussion

In this study, we simulate the temporal and vertical evolution of temperature and salinity in 19 individual mesocosms (PeECE III and KOSMOS 2013), by employing a one-dimensional water column model for mesocosms. The model has two key features that allow for realistic simulations of observed temperature and salinity patterns: (1) a density-based eddy diffusivity parametrization with one free parameter for optimization and (2) a surface correction of heat and freshwater fluxes, based on observed changes in total heat content and integrated water column salinity (see Sect. 2 for a detailed description). To minimize the model–data misfit, we optimize a diffusivity parametrization with respect to temperature and salinity. In this publication, we introduce our model, provide our best estimates of eddy diffusivities for each of the 19 mesocosms (see Supplement), and demonstrate that these estimates lead to simulated patterns of temperature and salinity that closely resemble the observed patterns (see Sect. 3 and Supplement).

4.1 Observed and simulated profiles of temperature and salinity

For both mesocosm experiments that we investigated in this study, PeECE III and KOSMOS 2013, our results show that vertical mixing patterns are very heterogeneous in time and space. For the PeECE III experiment, our simulations support the assumption of previous studies (Riebesell et al., 2007; Schulz et al., 2008) that water at the bottom of the PeECE III mesocosms was largely separated from the surface mixed layer during the whole experiment. Within the pycnocline that separates surface and bottom water, the simulated eddy diffusivities are very low, suggesting that tracer exchange between surface and bottom water was mostly inhibited. However, whether tracer diffusion through the pycnocline was still large enough to affect surface biogeochemistry can only be detected by a model that includes biogeochemical tracers explicitly. The PeECE III mesocosms show a distinct pattern for each optimization, due to the sensitive response of the artificial halocline to even slight changes in vertical diffusivities. In the KOSMOS 2013 experiment, the water column was fully mixed in the beginning and subsequently increasingly divided by thermal stratification. Our model is able to capture the very different mixing regimes, and it reproduces the observed temperature and salinity profiles accurately. In seven out of the 19 simulated mesocosms, either the temperature optimization or the salinity optimization results in the same diffusivity estimates as the optimization that depends on temperature and salinity simultaneously; i.e., in these cases the optimal diffusivities could have been estimated by taking into account only temperature or only salinity. However, since this is evident only in hindsight, after the different optimizations were conducted, we recommend always estimating diffusivities with the simultaneous temperature and salinity optimization, which provides the optimal diffusivities in any case.

4.2 Constant product of dissipation rate and mixing efficiency

For the calculation of eddy diffusivities, we keep the parameter c, the product of dissipation rate and mixing efficiency in the eddy diffusivity parametrization (Eq. 2), constant over time and space. All of our eddy diffusivity estimates are within the range of observed values found in the literature (e.g., Waterhouse et al., 2014; Whalen et al., 2012; Rovelli et al., 2016; Fer, 2009). The overall range of eddy diffusivities within a mesocosm is determined by the parameter c, which is optimized for each mesocosm independently. Thereby, the model can account for significant external influences, like local average wind speed or wave intensity. As we demonstrated in Sect. 3, eddy diffusivities in mesocosms of the KOSMOS 2013 experiment are much higher than those in PeECE III mesocosms, potentially due to differences in local wind speed, geographic position, and surrounding currents. We do not recommend using a constant c if significant seasonal changes in wind speed occur during the experiment. Both experiment sites we investigated here, Bergen and Gullmar Fjord, are characterized by a stormy winter season and a calmer summer season. Both mesocosm experiments were conducted within the calmer season.

4.3 Deriving diffusivities from observations

We derive eddy diffusivities from observed temperature and salinity profiles during the whole simulation period. Alternatively, it would be possible to only initialize the model with observed profiles and subsequently derive diffusivities from simulated temperature and salinity. However, our approach has the advantage that an accumulation of errors can largely be avoided. If diffusivities were derived from simulated temperature and salinity, even a small error in simulated values would lead to a progression of rapidly increasing errors. Furthermore, using observed temperature and salinity data also has the advantage that the influence of temporary storm events can be taken into account. A storm tends to increase mixing and decrease density gradients. Consequently, our model calculates higher diffusivities, since diffusivities and density gradients are inversely proportional (Eq. 2). This means, even though there is no direct dependency of diffusivities on wind speed, a temporary strong mixing event is automatically induced in the model by using observed temperature and salinity data.

4.4 Assumption of uncorrelated temperature and salinity

In the cost function that we use for the diffusivity optimizations, we assume that there is no correlation between salinity and temperature, because temperature and salinity can be positively or negatively correlated, depending on the environmental conditions. A strong increase in temperature at the surface would likely enhance evaporation, thus leading to an increase in salinity (positive correlation), a pattern that has been observed in the KOSMOS 2013 experiment. Conversely, in the case of thermohaline stratification, as present during the PeECE III experiment, warm water with a relatively low salinity lies on top of colder, more saline water – mixing of the two water masses would decrease surface temperature and increase surface salinity (negative correlation), and vice versa for water below the pycnocline. Since a combination of these processes could take place any time during the simulation period, temperature and salinity changes are treated as independent in the optimizations.

4.5 Limitations of the model

One important assumption in our model is that all changes in integrated water column heat and salinity are caused by heat and freshwater fluxes at the surface. The model does not account for heat fluxes through the mesocosm walls. In reality, heat transfer through the mesocosm walls can be significant when surrounding waters are influenced, e.g., by changing currents that transport colder water to the experiment site, as has been observed during the KOSMOS 2013 experiment (Bach et al., 2016). Presumably, our model overestimates heat in the lower part of the mesocosms because heat exchange with surrounding waters is unaccounted for. In the model heat can only exit the system after it is transported to the surface by turbulent diffusion. For freshwater fluxes, i.e., salinity changes, this limitation is not critical because even in reality there is no mass flow through mesocosm walls, as long as they are intact. In both experiments, temporary holes were detected but quickly fixed. Whether associated intrusion of surrounding waters had an impact on salinity inside the mesocosms is unknown. In any case, the assumption that heat (with the exception of penetrative solar heating) and freshwater fluxes occur only at the surface is a potential source for model errors. Another limitation of our model is the lack of nightly convection. Since temperatures were measured only at daytime, it is not possible to account for nightly cooling and related nightly convection. A continuous monitoring of temperature and salinity on site would allow for more exact diffusivity estimates. For future mesocosm experiments, we recommend putting data loggers into each mesocosm that record temperature and salinity continuously throughout the experiment. Thereby, nightly convection could be taken into account in a mesocosm mixing model and eddy diffusivity estimates could be further improved. Other processes that are not accounted for include reduction of solar radiation by wall shadowing and the direct effect of wind speed on turbulent mixing. Despite these limitations our simulated temperature and salinity profiles match the observations very well, suggesting that daily observations of temperature, salinity, and solar radiation are sufficient to simulate the main characteristics of vertical mixing patterns in mesocosms. For the purpose of biogeochemical modeling, the main characteristics like mixed-layer depth or pycnocline thickness are most important. Thus, the model we present here can be the basis for future biogeochemical mesocosm model studies.

4.6 Comparison with standard mixing assumptions

Most analyses of mesocosm data assume that there has been either no mixing or full mixing of the water column. To demonstrate that our model is significantly better than these simple assumptions, we compare the cost of our best estimate, J_TS(c_TS), with the cost of a simulation without any mixing other than convective mixing, J_TS(no mixing), and also with the cost of a simulation with full mixing, J_TS(full mixing).

In the case of the PeECE III mesocosms, when there is no mixing the cost J_TS(no mixing) ranges between 1.30 and 3.17, and the average cost of all mesocosms is 25 times higher than the cost of our temperature and salinity optimization J_TS(c_TS). For full mixing, the cost J_TS(full mixing) ranges between 0.53 and 0.55 and is on average 6 times higher than J_TS(c_TS).

In the case of the KOSMOS 2013 mesocosms, the cost J_TS(no mixing) ranges between 13.57 and 32.72 and is on average 302 times higher than J_TS(c_TS). For full mixing, the cost J_TS(full mixing) ranges between 0.07 and 0.29 and is on average 47 % higher than J_TS(c_TS).

Thus, we can confidently state that our mixing model is able to estimate much more realistic mixing conditions than simple assumptions of no mixing or full mixing and can reproduce observed patterns of temperature and salinity very well.

5 Code and data availability

The model code and the best estimates of the eddy diffusivities for all mesocosms of the PeECE III and KOSMOS 2013 experiment are available at https://doi.org/10.1594/PANGAEA.905311 (Mathesius et al., 2019).

6 Conclusions

This study sets out to simulate the physical processes in the mesocosm experiments PeECE III and KOSMOS 2013 (Schulz et al., 2008; Bach et al., 2016), with a focus on vertical mixing. We present here a one-dimensional mesocosm mixing model that is capable to reproduce the observed temporal evolution of temperature and salinity profiles of both mesocosm experiments despite very different environmental conditions. The density-based diffusivity parametrization that is used in our model yields plausible patterns of vertical diffusivities that are very heterogeneous over time and space. Optimal diffusivity estimates are achieved when the cost function depends on both temperature and salinity and the misfit between model results and observations is lowest. We provide here our optimal diffusivity estimates for the PeECE III and KOSMOS 2013 experiment (see Supplement). The development of the one-dimensional mixing model we present here can be regarded as an important foundation for future mesocosm model studies. Vertical mixing plays a significant role in mesocosms and has to be taken into account to accurately simulate observed temperature and salinity patterns, as we demonstrate in this study. Since our estimated eddy diffusivities can reproduce both temperature and salinity profiles well, we assume that these diffusivities can be applied to biogeochemical tracers as well. Including both vertical mixing and biogeochemistry in future mesocosm model studies could lead to a deeper understanding of biogeochemical processes in the pelagic zone, e.g., changes in carbon export due to ocean acidification as well as changes in primary production due to varying light exposure and nutrient availability. We encourage other modellers to work with the diffusivity estimates we provide here or to use our model to calculate diffusivities for any other mesocosm experiment, so that further insights can be gathered from observational data of mesocosm experiments. For future mesocosm experiments, we recommend measuring temperature and salinity profiles continuously with fixed data loggers within the mesocosms, so that the accuracy of diffusivity reanalyses can be further improved.