The influence of winter convection on primary production: a parameterisation using a hydrostatic three-dimensional biogeochemical model

In the recent past observational and modelling studies have shown that the vertical displacement of water parcels, and therefore, phytoplankton particles in regions of deep-reaching convection plays a key role in late winter/early spring primary production. The underlying mechanism describes how convection cells capture living phytoplankton cells and recurrently expose them to sunlight. This study presents a parameterisation called `phytoconvection' which focuses on the influence of convection on primary production. This parameterisation was implemented into a three-dimensional physical-biogeochemical model and applied to the Northwestern European Continental Shelf and areas of the adjacent Northeast Atlantic. The simulation was compared to a `conventional' parameterisation with respect to its influence on phytoplankton concentrations during the annual cycle and its effect on the carbon cycle. The simulation using the new parameterisation showed good agreement with observation data recorded during winter, whereas the reference simulation did not capture the observed phytoplankton concentrations. The new parameterisation had a strong influence on the carbon export through the sinking of particulate organic carbon. The carbon export during late winter/early spring significantly exceeded the export of the reference run. Furthermore, a non-hydrostatic convection model was used to evaluate the major assumption of the presented parameterisation which implies the matching of the mixed layer depth with the convective mixing depth. The applied mixed layer depth criterion principally overestimates the actual convective mixing depth. However, the results showed that this assumption is reasonable during late winter, while indicating a mismatch during spring.


Introduction
Photosynthesis in the ocean has been estimated to contribute almost 50% to the total global net primary production (Field et al., 1998) and total carbon uptake (Sabine et al., 2004). This large stake of marine primary production shows the great importance of understanding the physical, chemical and biological processes which influence marine primary production. Marine primary production is generally controlled by the availability of light and nutrients. In the ocean, light decreases exponentially with depth, limiting photosynthesis to the upper part of the ocean, the euphotic zone.
The 'critical depth model' (Sverdrup, 1953) describes the relation between the depth of the surface mixed layer and the capability of light-dependent phytoplankton net growth. It defines the compensation depth as the depth where the gain, or growth, and the loss in phytoplankton balance each other. Hence, there exists a critical depth where the vertically integrated growth is equal to the vertically integrated loss. Sverdrup (1953) concluded that net phytoplankton growth is only possible if the mixed layer depth (MLD) is less than the critical depth, thus, allowing for positive net growth.
During winter the MLD in the North Atlantic reaches depths of several hundred metres (e.g. McCartney and Talley, 1982) and therefore would not allow net phytoplankton growth according to the 'critical depth model', that is unless the loss terms are sufficiently low, and thus, result in a sufficiently deep critical depth (critical depth ≥ MLD). In this context, Behrenfeld (2010) proposed the so-called 'dilution-recoupling hypothesis', which based upon earlier work by Evans and Parslow (1985), arguing that the dilution of phytoplankton by mixed layer deepening reduces the grazing pressure on phytoplankton, and thus, lowers this loss term during winter. This hypothesis has been updated to the more general 'disturbance-recoupling hypothesis' (Behrenfeld et al., 2013;Behrenfeld and Boss, 2014) stating that during the annual cycle the predator-prey interaction is disrupted by environmental factors as, for example, convection, and subsequently recovers. This disruption then allows for net primary production due to low grazing pressure while concentrations still stay low due to convective mixing. Ocean convection -the buoyancy-driven sinking of surface waters due to surface cooling or an increase in surface salinity -is one of the key processes affecting the winter mixed layer deepening. Due to mass conservation the sinking of water parcels leads to a balancing upward motion, and thus, convection can be described as an orbital motion. Especially in the North Atlantic Ocean extensive convectively driven mixed layer deepening can be observed during winter. McCartney and Talley (1982) determined an average 'late winter' (January -April) MLD of more than 400 m over large parts of the subpolar North Atlantic using temperature and salinity profiles recorded during the 1950s and 1960s. Holliday et al. (2000) reported on convection in the Rockall Trough removing the seasonal thermocline and reaching down to depths of around 700 m.
In the meantime studies showed that convection plays a key role for primary production during winter (e.g. Backhaus et al., 1999;Wehde et al., 2001;Ward and Waniek, 2007). Backhaus et al. (1999) hypothesised a direct link between ocean convection and winter primary production called 'phytoconvection'. They suggested that the upward and downward motion within a convection cell causes phytoplankton particles to regularly re-enter the euphotic zone allowing them to grow, and thus, balancing their losses due to respiration, mortality, grazing and sinking. This relation was supported by modelling and observational studies (Wehde and Backhaus, 2000;Wehde et al., 2001;Backhaus et al., 2003). Ward and Waniek (2007) emphasise the role of convective upward motion which counteracts the sinking of phytoplankton, and thus, increases net growth.  to the Iceland Basin (e.g. Krauss, 1995) no stratification is visible in spring 2008. During late winter 2006 and 2007 the greatest MLD reached up to 600 m. The chlorophyll time series (B) is presented in relative units due to a lack of calibration data. However, in the chlorophyll the same patterns as in σ θ are shown underlining the importance of winter convective mixing for the phytoplankton community.
In large-scale ocean models, an adequate representation of the influence of convection on primary production is often missing. These models usually work on large horizontal and smaller vertical length scales, and therefore normally neglect the vertical acceleration. For many large-and mesoscale processes this does not cause any restrictions and the reduction of the vertical resolution saves computation time. However, when dealing with motions characterised by spatial aspect ratios of order 1, like convection, the hydrostatic approximation becomes inaccurate (Marshall et al., 1998). Furthermore, the coarse vertical resolution of these models cannot resolve the upward and downward motion of water parcels with high vertical velocities as it is known for ocean convection.
As a consequence of these two restrictions, models which do not a priori consider all variables to be homogeneously distributed within the mixed layer, may not represent the vertical exchange of water mass properties and phytoplankton in deep MLDs properly. This necessitates a parameterisation of primary production which is able to reproduce observed winter phytoplankton using these models (Holt et al., 2014).
Previous parameterisations of the influence of convection on primary production in Eulerian models lead to a reduction of primary production in deep MLDs (e.g. Lévy et al., 1998a,b). Furthermore, these studies focused on the onset of the bloom in spring after winter convection (Lévy et al., 1998b) or did not deal with mixed layers deeper than 200 m (Lévy et al., 1998a) while the study at hand focusses on primary production during winter when convection is present reaching depths of up to 500 m. Backhaus et al. (2003) suggested an approach for the parameterisation of phytoconvection taking into account the spatial aspect ratio of convective orbits. Janout (2003) adopted the idea and implemented it in a one-dimensional hydrostatic physical-biological model. The idea behind his approach was to compensate the lack of convective vertical displacement of phytoplankton by allowing primary production throughout the whole convective mixed layer. That was done by distributing the daily averaged surface solar radiation over the whole mixed layer during winter. To account for summer situations in which the 'conventional' approach following Sverdrup (1953) is more applicable, he switched between the conventional parameterisation and the phytoconvective approach depending on the MLD. This switching was performed abruptly by applying only the conventional parameterisation during periods with a MLD less than 75 m and applying the phytoconvective type for deeper mixed layers.
The present study incorporates the parameterisation of Janout (2003) and aims for the further development of this approach with two major objectives. First, to develop a smoother transition scheme between the conventional and the phytoconvective parameterisation of primary production. Second, to improve the applicability in a large-scale, three-dimensional, physicalbiogeochemical model applied to an area including regions of deep winter convection as well as shallow shelf seas.

The three-dimensional physical-biogeochemical model
For the three-dimensional (3D) simulations the ECOHAM4 model system (ECOlogical model, HAMburg, version 4;Pätsch and Kühn, 2008) was used. The model consists of two main modules: the hydrodynamic model HAMSOM (HAMburg Shelf Ocean Model; Backhaus, 1985) and the biogeochemical model ECOHAM. HAMSOM simulates the temperature and salinity distributions, the 3D advective flow field and the turbulent mixing. It is a baroclinic primitive equation model with a free surface and uses the hydrostatic and Boussinesq approximation. HAMSOM is defined on an Arakawa C-grid (Arakawa and Lamb, 1977) resolving the vertical in z-coordinates. The horizontal advective flow field is calculated using the component upstream scheme. A detailed description of HAMSOM is given by Backhaus and Hainbucher (1987) and Pohlmann (1991Pohlmann ( , 1996. The vertical turbulent mixing is parametrised by the exchange coefficient A V in depth z assuming stationarity and neglecting advection and diffusion of turbulent kinetic energy (Mellor and Yamada, 1974): Here, u and v represent the zonal and meridional velocity components, respectively. N is the Brunt-Väisälä-frequency and S m is the Schmidt number. H ml , representing the MLD, and C ml are the only terms which must be prescribed. C ml ≈ 0.05 is determined after Kochergin (1987). If unstable conditions occur (N 2 < 0), A V is set to the maximum vertical exchange coefficient A V max = 800 · 10 −4 m 2 s −1 to represent the strong vertical mixing due to convective processes which cannot be resolved by HAMSOM. The biogeochemical model ECOHAM describes the cycles of carbon (C), nitrogen (N), phosphorus (P), silicon (Si) and oxygen (O 2 ). It applies a partly variable C:N stoichiometry depending on the state variable and a simple description of benthic processes (Pätsch and Kühn, 2008). The model includes the following state variables: four nutrients (nitrate, ammonium, phosphate, silicate), two phytoplankton groups (diatoms and flagellates), two zooplankton groups (micro-and mesozooplankton), bacteria, two fractions of detritus (fast and slowly sinking), labile dissolved organic matter, semi-labile organic carbon, oxygen, calcite, dissolved inorganic carbon and total alkalinity. For a detailed description and a full list of the fluxes between the different state variables and the model parameters, see Lorkowski et al. (2012).

The parameterisation of 'phytoconvection'
In ECOHAM the phytoplankton production P B is controlled by an extended form of 'Liebig's law' (Liebig, 1840) considering light and nutrient limitation and accounting for the temperature dependence: The production rate P B depends on the maximum growth rate γ, the light L and the nutrients N . f T is a constant factor parameterising the temperature dependence of primary production. For our area of interest and the late winter situation it has been shown that light is the limiting factor during winter as the whole euphotic zone is enriched with nutrients due the upward mixing of nutrient-rich deep water driven by strong surface cooling, and thus, convection. Hence, in the following we will focus on the light dependence of phytoplankton production. The empirical relationship between irradiance and primary production found by Steele (1962) builds the basis for the conventional parameterisation of the light-dependent production P B at depth z: in which P B max depicts the maximum biomass-specific production rate (P B max = 1.1 d −1 for diatoms, P B max = 0.9 d −1 for flagellates) and I par is the depth-dependent photosynthetically active radiation (or PAR) and I opt is the dynamically calculated optimal light intensity. I par in depth z is calculated as: I par (z) = k par · I sw · exp(− (z)).
Here, k par = 0.43 (Pätsch and Kühn, 2008) represents the photosynthetically active fraction of the incoming short-wave radiation at the sea surface I sw . The depth-dependent attenuation coefficient (z) is: and includes the effect of light attenuation by water, planktonic self-shading and turbidity due to suspended particulate matter. k w is the locally varying attenuation coefficient of water after Jerlov (1976), k p = 0.03 and C p are the attenuation coefficient in m 2 mmol C −1 and the concentration in mmol C m −3 of phytoplankton, respectively. k s = 0.06 and C s are the attenuation coefficient m 2 g −1 and the concentration in g m −3 of suspended particulate matter, respectively. ECOHAM uses a variable I opt and its adaptation to the actual light conditions over time is described in Pätsch and Kühn (2008). The range of the optimum irradiance I opt is limited to: Starting from equation (3) the new parameterisation of light-dependent primary production under convective conditions is developed. Convection described as the upward and downward motion of water masses is here considered as an orbital system within the mixed layer. Turner (1979) and Kämpf and Backhaus (1998) studied the spatial and temporal scales of convective cells using observations and numerical simulations, respectively. The spatial aspect ratio (horizontal vs. vertical scale) has been reported to be on average about 2.5:1 (Wehde and Backhaus, 2000) and ranging between 2:1 and 3:1 (Kämpf and Backhaus, 1998).
The convection cell is considered having a rectangular geometry with the depth of the convective mixed layer H cml defining the vertical dimension (Backhaus et al., 2003). This represents a simplified approach to account for the spatial aspect ratio of convective orbits, since convective motion is highly turbulent and far from linear movement along a rectangular track. The conceptual view of such a convection cell is shown in Fig. 3. It is assumed that during winter the MLD is identical to H cml . Hence, the MLD determines the vertical extent of the convection cell. Using the spatial aspect ratio of 2.5:1 and the time of a complete orbit t orb can be calculated as (Backhaus et al., 2003;Janout, 2003): Here, v c is the velocity of the convective motion which is assumed to be constant with v c = 5 cm s −1 . Following Janout (2003), who adapted Backhaus et al. (2003) the time within the euphotic zone is defined as: with the euphotic depth H euph . H euph is set to the depth where the available light is still 1% of the surface radiation. H cml is defined as the last depth z in which the temperature difference criterion (9) is satisfied.
The value of ∆T = 0.4 K is within the range of literature values varying between ∆T = 0.1 K and ∆T = 1 K compared to the sea surface temperature (SST) (see e.g. Table 1 in Kara et al. (2000)).
Using equations (7), (8) and the MLD criterion (9) the ratio t exp /t orb can be calculated. This ratio implies that each phytoplankton particle within a convective cell has the same probability of residence within the euphotic zone which is the fundamental idea behind the parameterisation presented here. It is assumed that phytoplankton production still follows equation (3) and is conducted in the euphotic zone with the average production rate under the actual light conditions. Combining this assumption with the above ratio t exp /t orb , we estimate the 'phytoconvective production' P B pc : which is constant throughout the whole mixed layer. It follows that under convective conditions each plankton particle has the same probability to conduct primary production with the average light within the euphotic zone independent of its actual position within the mixed layer. The only term changing in equation (10) is the ratio t exp /t orb , since equation (10) is only applied if H cml > H euph and P B Steele and H euph are assumed to be identical for different H cml . Thus, even though t exp increases with increasing H cml , t exp /t orb decreases because of the faster increase of t orb . In other words, the frequency of convective orbits is inversely correlated with H cml and the same holds for the exposure to light (Backhaus et al., 2003), causing P B pc to be reduced with increasing H cml . Convection of several hundreds of metres depth only occurs during winter, making it necessary to distinguish between convective and non-convective periods. The transition between convective and non-convective regimes is controlled by H cml which is used to build the weighting function (11). This weighting function controls the influence of the conventional parameterisation P B Steele (equation (3)) and the phytoconvective parameterisation P B pc (equation (10)). Furthermore, there is a discontinuity between P B Steele and P B pc in the case of H cml = H euph due to the additional distance 2.5 · H cml travelled at the surface (see Fig. 3). This discontinuity is eluded by the use of the weighting function.
This weighting function f p is set to 0, if H cml is less than or equal to the euphotic depth H euph . Accordingly,  (10).
The weighting factor f p is then applied for the calculation of the actual light-dependent phytoplankton production in depth z: During summer, when the mixed layer is shallower than the euphotic zone, the phytoconvective parameterisation (10) is not taken into account. Conversely, during winter, when the MLD reaches several hundreds of metres the conventional parameterisation (3) is not affecting the phytoplankton productivity. The light-dependent production rate P B total according to equation (12) is consequently used for the calculation of primary production within the model.

Model setup and data
The model was set up on the region of the Northwestern European Continental Shelf (NECS) and parts of the adjacent Northeast Atlantic (Lorkowski et al., 2012). The resolved model region is shown in Fig. 4 We focused our analysis of the new parameterisation on a station west of Ireland in the eastern Rockall Trough (55 • 3' 24" N, 10 • 14' W, see marked position in Fig. 4). This station was selected because winter convection reaching depths to more than 500 m has regularly been reported in the Rockall Trough (Holliday et al., 2000). Meincke (1986) reported MLD of up to 1000 m during severe winters. Furthermore, it is one of the few stations where chlorophyll measurements are available for the winter season. The observation data was provided by the British Oceanographic Data Centre (BODC) and was measured on February 27th, 1996, during the cruise CH125B of RRS Challenger (Hill, 1996). The hydrodynamical model HAMSOM was initialised with a monthlyaveraged climatology based on the World Ocean Atlas (WOA; Conkright et al., 2002). Salinity was treated as a semi-prognostic variable and adjusted to the climatology with a time constant of 14 days to guarantee reasonable salinity distributions. At the open boundaries temperature and salinity are prescribed during inflow situations with a time constant of 7 days. Additionally, the surface elevation according to the M2 tide was prescribed at the open boundaries. The meteorological forcing including air temperature, cloud coverage, relative humidity, wind speed and direction was calculated from NCEP/NCAR reanalysis data (Kalnay et al., 1996) and was provided as 6-hourly values. Starting from the initial data, the HAMSOM model applied a 10 day spin-up during which the air pressure was constant and the wind speed was zero. The results of the hydrodynamical simulation for 1996 were calculated as averages over two tidal cycles representing daily values. This simulation output served as the basis for the biogeochemical simulation with ECOHAM.
To demonstrate the effect of the new parameterisation, we compared a simulation run using the new parameterisation (hereafter 'phytoconvection run') to a simulation using the conventional parameterisation (hereafter 'standard run'), using equation (12) and equation (3), respectively. For both simulations, the initialisation of the biogeochemical state variables was done using a dataset produced by a model run using only the conventional parameterisation following equation (3). All other settings were according to Lorkowski et al. (2012).

The convection model
A non-hydrostatic convection model was used to evaluate the assumption that the MLD is persistently mixed by convection, and thus, constitutes a valid indicator for the vertical extent of the convective cell. It was furthermore used to test the validity of the ratio of t exp /t orb . The applied convection model uses the Boussinesq equations for an incompressible fluid in a 2,5-dimensional ocean slice with cycling boundary conditions, thus, allowing to include the rotational effects of the earth. For the turbulent eddy viscosity the turbulence closure scheme by Kochergin (1987) is used. The model calculates the density (temperature and salinity), the hydrodynamic flow fields and the non-hydrostatic pressure over space and time. It neglects the influence of wind stresses but allows latent and sensible heat fluxes due to fluctuations in the wind speed. Thus, a free convective boundary layer is assumed. For a detailed description of the model, including the equations, the reader is referred to Kämpf and Backhaus (1998) and to Wehde and Backhaus (2000).
The model was set up on a domain of 1000 m depth and 250 m width with a horizontal and vertical grid size of 5 m applying a time step of 10 seconds. We conducted two simulations starting at March 17th and March 29th to capture the transition period between winter convective mixing and the decline of the mixed layer, respectively. Each simulation ran over a period of 14 days with one additional day as spin-up. The initial profiles of temperature and salinity for the first simulation were vertically interpolated from the simulated profiles of the 3D simulation at the station in the Rockall Trough (see Fig. 4). The second simulation was initialised with the results of the first one. The same meteorological forcing as for the 3D simulation for the respective station was used, provided as 3-hourly interpolated data. To test the validity of the use of the ratio of t exp /t orb 200 Lagrangian tracers were randomly distributed within the mixed layer at the beginning of the first simulation. They record the actual light throughout the simulation period, thus, allowing the calculation of t exp and t orb for each tracer.

Results and Discussion
3.1. The three-dimensional physical-biogeochemical model We will focus our discussion on the biogeochemical simulations as the physical simulation only built the basis for testing the parameterisation of phytoconvection. For this purpose, only a brief presentation of the simulated temperature and the derived MLD is given. All Hovmöller diagrams, time series and the comparison to the observations refer to the station in the Rockall Trough marked in Fig. 4.

Hydrodynamic simulation
The temporal development of the simulated temperature at the station in the Rockall Trough (see Fig. 4) and the corresponding MLD (dashed lines) are shown in Fig. 5. The MLD resulting from equation (9) (black line) is compared to the MLD determined using a temperature difference (dark grey line) and density difference criterion (light grey line) after Kara et al. (2000) based on a ∆T of 0.8 K. The simulation starts from the climatological dataset with mixed layer temperatures of above 10 • C which steadily decrease until mid-March due to ongoing surface cooling and deep mixing. During this period the MLD determined after equation (9) is at maximum 500 m which is in the same order as the climatological MLD reported for the northern Rockall Trough at about 57 • N to 58 • N (Holliday et al., 2000), which is at maximum about 700 m deep. In April the mixed layer declines and a shallow seasonal surface mixed layer develops reaching maximum temperatures of above 13.5 • C in August. Starting in October, the MLD deepens again reaching depths of up to 200 m at the end of the year. Throughout the whole simulation period the two different MLDs determined after Kara et al. (2000) show greater values than the applied criterion. From January to March the two criteria yield MLDs 100 m to 300 m deeper than the MLD criterion used in our study, which contradict the temperatures visually decreasing already in 300 m to 400 m depth. In mid-April the strong decrease in the applied MLD coincides with the beginning of near-surface thermal stratification whereas the two MLDs after Kara et al. (2000) stay deep until late April showing. This supports the use of the applied MLD criterion.  (9), the dark grey and light grey lines refer to the MLD determined with a temperature difference and, respectively, density difference criterion based on ∆T = 0.8 K (Kara et al., 2000). The monthly averaged MLD after (9) for March and August within the whole model area is shown in Fig. 6. The average MLD of 400 m to 600 m in March in the Rockall Trough region are in good agreement with Holliday et al. (2000). In most parts of the shelf the MLD represents the bottom topography due to strong winter mixing reaching down to the bottom. In August most parts of the model area are stratified with MLD less than 20 m representing a small underestimation (Pätsch and Kühn, 2008). Only parts of the English Channel, the Irish Sea and the southern North Sea show a MLD reaching the bottom. This is in good agreement with Pingree and Griffiths (1978) who reported on perennial well-mixed waters in these regions due to strong tidal mixing. Fig. 7 shows the monthly and vertically averaged specific light limitation function (LLF) within the mixed layer for the standard run (dark grey) and the phytoconvection run (light grey) at the station west of Ireland (see Fig. 4). The LLF is defined as the light-dependent production rate (equation (12)) normalised by the maximum production rate P B max , and therefore, is a direct measure for the phytoplankton production.  Figure 7: Time series of the monthly and vertically averaged specific light limitation function within the mixed layer for the standard run (dark grey) and the phytoconvection run (light grey). The specific light limitation function is defined as the light-dependent production rate (see equation (12)) normalised by the maximum production rate P B max .

Biogeochemical simulations
The two simulations generally follow the seasonal cycle of the solar radiation showing low values during winter and high values during summer. Because of the maximum influence of the phytoconvective parameterisation, the phytoconvection run shows significantly higher values than the standard run during winter due to the higher values of the LLF in the deeper part of the mixed layer. From May to October the two simulations show similar LLFs, due to the shallow seasonal mixed layer causing the influence of phytoconvection to be zero most times. Small differences between the two simulations during this period are caused by the interaction of the conventional parameterisation and the phytoconvective parameterisation during events of wind-induced mixed layer deepening. At the end of the year the LLFs are again diverging due to the deepening of the mixed layer, and thus, the increasing influence of phytoconvection. Fig. 8 shows the Hovmöller diagrams of chlorophyll-a (chl-a) simulated by the standard run (A) and the phytoconvection run (B) at the station west of Ireland. The chl-a concentrations were derived from the simulated phytoplankton carbon applying a constant chl-a:C mass ratio of 1:50. The phytoplankton carbon is the sum of carbon of the two simulated phytoplankton groups diatoms and flagellates. The two simulations start from the same initial conditions with highest concentrations of above 0.1 mg chl-a m −3 in the upper 100 m and decreasing concentrations in greater depths.
In the first two weeks, the two simulations show decreasing concentrations in the upper 100 m and, due to downward mixing of the higher surface concentrations, increasing concentrations in the layers in between 100 m and the MLD. Thereafter, the chl-a of the two simulations diverge. The standard run shows a steady decrease until late March throughout the whole water column while the phytoconvection run is characterised by steadily increasing concentrations throughout the mixed layer from early February until mid-April. This development is controlled by the different LLFs Fig. 7. The concentrations in the standard run drop below 0.05 mg chl-a m −3 in the upper 100 m and values less than 0.035 mg chl-a m −3 in the deeper layers, while the phytoconvection run shows concentrations of above 0.05 mg chl-a m −3 throughout the whole mixed layer and concentrations higher than 0.1 mg chl-a m −3 in the upper 200 m to 300 m. Maximum concentrations in the deep part of the winter mixed layer are reached in early April when the MLD declines and a seasonal thermocline develops with values of above 0.35 mg chl-a m −3 in the phytoconvection run. At the same time the standard run shows the onset of a strong surface spring bloom indicated by the increase in the chl-a concentrations from less than 0.05 mg chl-a m −3 to above 1 mg chl-a m −3 within about two weeks. A comparably strong surface bloom is not simulated in the phytoconvection run, because the MLD is still 200 m to 400 m deep, and hence, the phytoconvective parameterisation is fully taken into account.
From late April to November the two simulations show a very similar behaviour within the mixed layer as expected from the similarity of the LLFs during this period (see Fig. 7). Caused by the increasing influence of the phytoconvective parameterisation in December, the phytoconvection run shows slightly higher concentrations than the standard run in the deeper layers (50 m to 300 m). Fig. 8 showed that the phytoconvection run simulates significantly different winter chl-a concentrations than the standard run. In order to evaluate the simulated chl-a profiles, the standard run (dark grey, dash-dotted) and the phytoconvection run (dark grey, dashed) are compared to two measured chl-a profiles (solid) recorded on February 27th, 1996 (Hill, 1996) (Fig. 9).
The horizontal black lines mark the MLD determined using equation (9) from the simulated (dashed) and measured (solid) temperature profiles, respectively. The two observed profiles show a distinct structure with increased chl-a concentrations ranging between 0.09 mg chl-a m −3 and 0.14 mg chl-a m −3 in the upper 500 m to 600 m. Thereunder, concentrations strongly decrease followed by concentrations of about 0.05 mg chl-a m −3 below the MLD. Chla concentrations in the upper hundreds of metres are in the same order of magnitude as measurements made in different regions of the North Atlantic (Backhaus et al., 2003). The standard run shows maximum concentrations of about 0.05 mg chla m −3 at the surface and steadily decreasing concentrations with increasing depth. In contrast, the phytoconvection run shows increased concentrations of about 0.11 mg chl-a m −3 to 0.12 mg chl-a m −3 in the upper 200 m to 300 m followed by a stronger decrease and concentrations less than 0.02 mg chl-a m −3 below 700 m. In the upper 200 m to 300 m these results are comparable with the observations, even though the vertical variability of the measurements is not reproduced. The standard run shows significantly lower chl-a concentrations throughout the whole water column and does not reflect the distinct vertical structure of the observations. The vertical structure of the phytoconvection run fits the observations much better with its higher concentrations in the upper hundreds of metres and the subsequent stronger gradient between 300 m and 700 m with the maximum gradient around the MLD. Nevertheless, the simulated chl-a is not homogeneously distributed throughout the whole mixed layer which may be due to the vertical mixing of the phytoplankton to deeper layers.
In the observations the determined MLD is deeper than the depth of the maximum gradient in the chl-a which indicates that the applied MLD criterion may overestimate the actual MLD. In the simulation the maximum gradient corresponds to the simulated MLD. However, it remains constrained by the vertical resolution of the model grid which is only 100 m between 200 m and 600 m, and therefore is not able to reflect the sharp gradient of the observations.
This comparison to observations is rather basic due to the low data availability and gives just a first insight in the capabilities of the new parameterisation. Hence, a more compelling validation is required for a complete evaluation of the presented parameterisation, which we leave to future work.
The comparison of the chl-a concentrations with those in Lévy et al. (1998a) in February (the month with maximum MLD in their study) shows that our parameterisation yields concentrations in the same order of magnitude (about 0.1 mg chl-a m −3 ) while chl-a concentrations at the beginning of the period simulated by Lévy et al. (1998a) are about twice as high compared to our simulation. Thus, the parameterisation by Lévy et al. (1998a,b) applied on our initial conditions wouldmost likely result in chl-a concentrations underestimating the observations by about factor 2. Furthermore, the approach presented here is reversed to the approach by Lévy et al. (1998a,b) regarding the effect of convection on primary production. In our parameterisation primary production is increased by convection while Lévy et al. (1998a,b) included convective mixing as a limiting factor for primary production. The parameterisation by Lévy et al. (1998a,b) applies an additional factor γ m ranging between 0.1 and 1 representing maximum limitation of growth due to vertical mixing in the case of H cml > 2 · H euph and no convection-induced limitation for H cml < H euph , respectively. However, this parameterisation has not been applied to mixed layers deeper than 200 m and Lévy et al. (1998a) do not deal with cases of H cml > 2 · H euph . Thus, the win-ter chl-a concentrations simulated by Lévy et al. (1998a) would be reduced further in the case of a winter MLD of similar depth as in our simulations. Consequently, the parameterisation by Lévy et al. (1998a,b) is unlikely to reproduce the observed winter chl-a concentrations. We use the vertically averaged LLF which is subsequently distributed over the MLD whereas Lévy et al. (1998a,b) used the vertically averaged PAR. The advantage of using the vertically averaged LLF is, that it allows for photoinhibition which might occur in the surface layer when solar radiation increases in early spring. In contrast, using the average PAR reduces the PAR in the upper layers, and thus, can prevent photoinhibition resulting in an overestimation of growth.
The analysis of the phytoconvection run with respect to nutrients (specifically nitrate, not presented here) showed that even though the increased primary production during winter affects the vertical nutrient distribution, the levels stay well above any limiting thresholds. For the zooplankton grazing (also not presented) there are some differences between the standard run and the phytoconvection run, especially in the period from February to April. The monthly and vertically integrated grazing rates in the phytoconvection run exceed those in the standard run by above factor 8 during March. However, during the bloom initiation in mid-April the daily zooplankton grazing (not shown here) in the surface layer is very similar in both simulations (maximum deviations of about 13% with even higher rates in the standard run), indicating only a small impact on the phytoplankton development.
In ECOHAM, zooplankton grazing is directly proportional to the phytoplankton concentration which explains the stronger grazing in the phytoconvection run during winter. The model implicitly accounts for the decoupling of zooplankton and phytoplankton as zooplankton grazing declines with decreasing phytoplankton concentrations, and consequently zooplankton concentrations decrease. These dynamics are in agreement with the dilution phase of the 'Dilution-Recoupling-Hypothesis' (Behrenfeld, 2010), representing one possible approach for how winter phytoplankton concentrations can be reproduced in models. The 'recoupling phase' (Evans and Parslow, 1985;Behrenfeld, 2010) however is not captured by the model, since zooplankton, forced by the vertical exchange coefficient A V (see equation (1)), is mixed within the water column like other state variables, and thus, is unable to retain its position within the water column. To analyse the functioning of the developed weighting function (11) Fig. 10 shows the vertically integrated, monthly averaged chl-a concentrations of the upper 500 m for March (A and B) and August (C and D). In March the standard run (A) shows significantly lower concentrations than the phytoconvection run (B) in the areas south and west of Ireland due to the MLD being deeper than the reference depth H ref = 100 m (see Fig. 6). The highest concentrations in these regions exceed 600 mg chl-a m −2 in the phytoconvection run while the concentrations in the standard run stay below 150 mg chl-a m −2 in the same area.
In the deep areas north of the shelf the concentrations in both simulations are similarly low despite the deep mixed layers. Here, phytoplankton concentrations are very low at the beginning of the simulation (about factor 20 less than at the analysed station) and stay low throughout the winter. South of Norway both simulations show increased chl-a concentrations, but with a spatial mismatch. Due to the deep mixed layer the phytoconvection run produces the highest concentrations in the deepest area of the Skagerrak. In contrast, the standard run shows the highest concentrations directly at the coast due to the onset of the spring bloom. In the shallower shelf regions, e.g. the southeastern North Sea the two simulations are in good agreement with each other demonstrating that the weighting function allows the application of the new parameterisation on the shelf.
In August (C and D) the two simulations show the same patterns and concentrations in most parts of the shelf and the areas off the shelf. The shallow seasonal mixed layer in these areas (see Fig. 6) causes phytoconvection in the phytoconvection run (D) to switch off. Only in the English Channel and the Irish Sea the phytoconvection run shows significantly higher concentrations than the standard run (C). In these regions the bottom depth is around 100 m and the vertical mixing is strong throughout the whole year due to tidal mixing and the interaction of horizontal currents with the bottom topography. This prevents the development of a persistent seasonal mixed layer, i.e. the MLD deepens frequently causing the influence of phytoconvection to increase. Consequently, the chl-a concentrations in the deeper layers are significantly higher than in the standard run resulting in higher vertically integrated concentrations. These regions are also known for the regular occurrence of tidal fronts and increased phytoplankton production during summer due to the availability of light and nutrients . This indicates reasonable chl-a concentrations simulated in the phytoconvection run in these regions.
The monthly and vertically integrated (0 m to 500 m) net primary production (NPP) in the whole model region shows similar patterns as the chl-a ( Fig. 10) (not shown). In regions where the influence of phytoconvection is strongest (southwest of Ireland) the NPP in the phytoconvection run exceeds the NPP in the standard run by about factor 8 to 10 during March. Maximum values of above 18 g C m −2 month −1 are reached in these areas in the phytoconvection run. In August, the phytoconvection run reaches maximum NPP of above 30 g C m −2 month −1 in the English Channel while a maximum NPP of about 16 g C m −2 month −1 is reached south of the Doggerbank in the standard run.
To illustrate the effect of the new parameterisation on the carbon cycle Fig. 11 shows the monthly integrated air-sea carbon flux (ASCF, A) and the carbon export production (CEP, B). Negative values in the ASCF imply outgassing of CO 2 . The CEP is defined as the export of fast-sinking detritus below 500 m depth as this is the maximum MLD at the specific station. For the ASCF the two simulations are only slightly different from each other throughout the whole seasonal cycle. During winter the two simulations show maximum outgassing of CO 2 , or the most negative ASCF due to strong mixing, and hence, the transport of carbon-enriched deep water to the surface. During this time the phytoconvection run (light grey) shows slightly higher values induced by the higher phytoplankton biomass, and therefore, increased primary production at the surface compared to the standard run (dark grey).  Figure 11: Time series of monthly integrated (A) air-sea carbon flux and (B) carbon export production for the standard run (dark grey) and the phytoconvection run (light grey). The negative values in the air-sea flux mean outgassing of CO 2 . The carbon export production is defined as the export of fast-sinking detritus below 500 m depth referring to the maximum MLD at the analysed station.
During the spring bloom formation, the ASCF increases significantly due to the increasing near-surface primary production. As near-surface primary production during April and May is stronger in the standard run than in the phytoconvection run, the relation between the two simulations shifts and the ASCF becomes higher in the standard run. During summer the ASCF stays on a relatively high level due to the ongoing primary production. The slightly lower values in the phytoconvection run are caused by wind-induced mixed layer deepening, and thus, the influence of phytoconvection which lowers the primary production in the surface layer. From October to December outgassing again intensifies due to increased mixing which brings carbon-enriched water to the surface, and the reduced primary production. The standard run shows only slightly higher values. This is caused by the increased importance of the phytoconvection leading to lower near-surface production rates, while the near-surface phytoplankton biomass remains sim-ilar in the two simulations after summer (see Fig. 8).
While the ASCF is very similar for the two simulations throughout the year, the CEP shows significant differences during the annual cycle. In the standard run the CEP is steadily decreasing from January to April when the spring bloom is fully developed due to the decreasing phytoplankton biomass. In contrast, the phytoconvection run shows a steady increase in the CEP from February to April with maximum amounts being 6 to 9 times higher than in the standard run. This is caused by the significantly higher biomass of phytoplankton and zooplankton in the water column from the surface to the MLD, and consequently, higher amounts of fast-sinking detritus at great depths. The peak in April with a CEP of about 18.5 mmol C m −2 is related to the strong decline of the mixed layer leaving a large portion of the phytoplankton biomass below the mixed layer which is than transferred to detritus sinking to the deep ocean. During summer the two simulations are again converging and show a similar CEP until October. Thereafter, the CEP becomes again slightly higher in the phytoconvection run due to the increasing influence of phytoconvection, and hence, the slightly higher concentrations of phytoplankton in greater depths (see Fig. 8).
As the spatial aspect ratio of the convection cell (see Fig. 3) is a key factor for the parameterisation of phytoconvection we compared the chl-a simulated by the use of three different aspect ratios. Fig. 9 shows the effect of the applied spatial aspect ratios of the convection cell on the phytoplankton. The previously discussed phytoconvection run (dark grey, dashed) applying an aspect ratio of 2.5:1 is compared to simulations applying aspect ratios of 2:1 (light grey, dash-dotted) and 3:1 (light grey, dash-dotted), respectively. The presented ratios cover the range reported by Kämpf and Backhaus (1998). Simulated chl-a concentrations increase with rising aspect ratios as the time within the euphotic zone t exp increases (see equation (8)). Furthermore, it can be seen that for the first observed profile (left panel) the 2:1 ratio produces the best fit regarding the concentrations in the upper 300 m, while for the second observation (right panel) the 3:1 ratio fits best to the observations. This suggests, that on average the aspect ratio of 2.5:1 offers a reasonable fit.
The deviation in the simulated chl-a of the 2:1 and 3:1 ratio simulation relative to the phytoconvection run for the purely phytoconvective period (January to March) was calculated for each time step and layer at the presented station as the Euclidean distance normalised by the average of the the two simulations. The maximum deviation for the simulation applying a ratio of 2:1 is 23.7% reached on March 30th. The 3:1 simulation reaches its maximum of 19.8% on March 26th. These deviations are reached when the chl-a concentrations start to significantly increase (compare Fig. 8) due to enhanced light availability. These deviations would make a more thorough validation of the results originating from different spatial aspect ratios desirable. However, this would require a more comprehensive dataset, which unfortunately was not available.
The reference depth H ref controlling the weighting of the standard parameterisation and the phytoconvective parameterisation (see equations (11) and (12)

The convection model
In order to validate our assumption that the MLD (after equation (9)) is generally a good indicator for the convective mixing depth H cml , we applied the convection model described in section 2.2 to two different periods: March 17th to March 31st and March 30th to April 13th. Fig. 12 shows the merged time series of simulated temperature (A) and turbulence as vertical mixing coefficient (B) from March 17th to April 12th. The colour scale for the turbulence is cut at its upper end at 8 cm 2 s −1 to better resolve the values in greater depth. The vertical dashed lines (March 30th) mark the date of assembly. The daily running average of the MLD determined from the simulated temperature using the MLD criterion (9)  The temperature shows a steady decrease within the upper 400 m from the beginning of the simulation until April 2nd (day 93). This is in good agreement with the convective turnover times of 5 hours to 8 hours indicated by increased turbulence of about 2 cm 2 s −1 down to depths of 300 m to 400 m during this period. The MLD in principle follows the temporal development of the turbulence, but the MLD is estimated about 150 m deeper than the maximum depth of convective mixing during the whole simulation period. While the applied MLD criterion overestimates the H cml , the MLD determined from the convection model and from the 3D simulation are in good agreement (see Fig. 5). The simulated turbulence also shows that the MLD is recurrently fully convected during winter, thus, providing a good indicator for the vertical extent of the convection cell.
In the period after April 2nd, convection is strongly reduced. These dynamics are also visible in the temperature increase near the surface. During the same time, the MLD shows only a slight decline demonstrating that the MLD and the H cml deviate stronger during early spring. This period of three days of no convective mixing is followed by another period of strong convection causing a deepening of the mixed layer by almost 100 m. The temperature within the whole mixed layer also significantly decreases during this time. Convection eventually drops on April 6th (day 97) followed by the formation of a shallow surface mixed layer. The onset of thermal stratification is delayed compared to the shutdown of convection by about two days. This shows that especially in the time directly after the shutdown of convection the MLD is not a good indicator for the H cml for the relevant dynamics to describe phytoplankton growth conditions. This is also in good agreement with Ferrari (2011) andFerrari et al. (2014) who stated that during periods of strong atmospheric forcing the MLD is a good proxy for H cml , whereas the MLD strongly deviates from H cml when the atmospheric forcing weakens in spring, thus, becoming a less strong driver for turbulent mixing. However, as our study focusses on primary production during winter when the MLD is frequently overturned by convection, we do not further discuss this problem here.
The ratio t exp /t orb is about 0.38 throughout the winter period (January -March) in the 3D phytoconvection run due to the spatial aspect ratio of 2.5:1 (horizontal to vertical) of the convection cell and an euphotic layer depth of 40 m to 50 m. To evaluate this value we calculated the median values of t exp and t orb for the 200 Lagrangian tracers added to the first simulation. The median was used rather than the average value to retrieve values less influenced by strong outliers occurring in the tracer results. We derived the values by considering tracers with an actual light value of above or equal 1% of the sur-face light being within the euphotic zone, thus, counting for t exp . Tracers with light values below this threshold were taken into account for the dark time t dark . t orb was calculated as the sum of the median values of t exp and t dark . The resulting median values for t exp and t orb are about 0.08 hours and about 0.42 hours, respectively, resulting in a ratio t exp /t orb of about 0.2. This shows that the ratio t exp /t orb during the 3D simulation is most likely overestimated causing an overestimation of the simulated phytoplankton. The short median duration of t exp and t orb suggests that particles in a convectively mixed layer are not always transported down to the bottom of the mixed layer, but spend a lot of time in the upper layers and frequently re-enter the euphotic zone. However, laboratory experiments showed that phytoplankton dealing with short light-dark cycles is less productive than that dealing with longer cycle due to photoacclimation (B. Walter, unpublished data). This is not taken into account in our parameterisation as the different light-dark cycles are not represented and would probably lower the simulated phytoplankton concentrations. However, all tracer trajectories (not presented) show that particles also spend longer continuous periods in different depths since convection is intermittent. Thus, it is likely that a certain proportion of phytoplankton stays sufficiently long in the euphotic zone and net primary production can take place.
There are some differences between 2D and 3D modelling of turbulent convection, a point worth noting. For instance, Moeng et al. (2004) showed that in the case of a free convective boundary layer as assumed by the model, the turbulent kinetic energy (TKE), surface friction velocity and velocity variances are sensitive to the subgrid-scale eddy viscosity and thermal diffusivity in 2D models. Especially the vertical velocity variance is significantly higher in 2D models than in 3D which affects the vertical mass fluxes (Moeng et al., 2004). Hence, a 3D simulation using the same boundary conditions would lead to different results regarding the TKE, and thus, to a different H cml and MLD, respectively. Fox-Kemper and Ferrari (2008) showed that at the mixed layer base 2D models tend to simulate a lower vertical diffusivity compared to 3D models which may affect the temporal development of H cml . Even though such model would most likely yield a quantitatively different result, this does not affect the validity of our assumption that the MLD equals the H cml during winter as we do not expect another qualitative result from a 3D model. With respect to the vertical displacement of particles by convection, the differences in the vertical velocities would affect the simulated duration of a full convective orbit, and thus, the ratio t exp /t orb .

Conclusion and perspectives
The presented parameterisation of phytoconvection (equation (10)) combined with the developed weighting function (11) constitutes a possible approach to reproduce observed winter phytoplankton concentrations (e.g. Backhaus et al., 2003) in a large-scale physical-biogeochemical model. In these models the vertical exchange is often not sufficient to account for convective mixing, and thus, to adequately simulate phytoplankton concentrations within the deep winter mixed layer. This parameterisation incorporates the hypothesis of Backhaus et al. (1999) by allowing primary production throughout the whole mixed layer, and thus, increasing winter phytoplankton productivity. The resulting parameterisation (equation (12)) is able to simulate winter phytoplankton concentrations in convective regions which are in good agreement to the available observations whereas the conventional parameterisation of primary production fails to reproduce the observations. The more realistic winter phytoplankton concentrations and the accordingly higher carbon export production in the phytoconvection run suggest that convection non-resolving, biogeochemical models underestimate the carbon export as they also underestimate the phytoplankton concentrations in these depths, and hence, the sinking of particulate organic carbon.
Beyond this, for the first time the parameterisation proposed by Backhaus et al. (2003) was tested and implemented into a large-scale 3D biogeochemical ocean model. Thus, it was possible to investigate the effect of this parameterisation on the simulated ecosystem, for example, nutrients and zooplankton. Even though only primary production is directly influenced by the parameterisation, nutrients and zooplankton are also affected. However, the increase in zooplankton grazing and the reduction of nutrients in the pre-bloom phase only have a minor effect on the phytoplankton bloom development. Thus, these results question the relative importance of changes in the density-dependent grazing pressure during winter as formulated in the 'disturbance-recoupling hypothesis' (Behrenfeld et al., 2013;Behrenfeld and Boss, 2014).
In addition, the regional model ECOHAM4 served as a testing environment for the application of the parameterisation of phytoconvection in a global ocean model. The promising results with respect to phytoplankton concentrations suggest this step to investigate the effect on the phytoplankton and carbon export on a global scale.
The parameterisation presented relies on a number of assumptions re-garding the characteristics of convection and the functioning of primary production. First, it is assumed that the MLD is frequently mixed by convection during winter, thus, providing a valid measure for the vertical extent of convection. This was confirmed by the results of the convection model, although the applied MLD criterion overestimates the actual convective mixing depth H cml . Second, the convection cell is assumed to have a rectangular shape with a horizontal to vertical aspect ratio of 2.5:1 (Backhaus et al., 2003) and particles moving with constant velocity along linear pathways. The median timescales of t exp and t orb resulting from the tracer experiment showed that this simplification leads to an overestimation of the ratio t exp /t orb . The short timescales also suggest that the actual motion within the convection cell is rather turbulent than linear which also effects the particle velocity. Horizontal current velocities in the ocean are usually larger than vertical velocities. Hence, assuming increased horizontal velocities would lead to a reduction of t exp and consequently lower primary production. The short light-dark cycles show that the simulated phytoplankton is biased by neglecting the effect of photoacclimation resulting in a further overestimation. Consequently, taking into account these aspects would most likely result in an overall lower primary production during winter. Since primary production depends on absolute values of PAR, primary production taking place in depths deeper than the 1%-threshold depth could partly balance this decrease in primary production as it would increase t exp . Though it is not likely to have a strong effect during mid-winter when radiation is generally low, it may have a beneficial effect during March when surface radiation increases, and thus, primary production may occur in greater depth. Variation in respiration is another aspect not accounted for in the applied model. It is known that phytoplankton respiration rates decrease under low light conditions (e.g. Falkowski et al., 1985). Thus, including this effect is expected to increase the simulated NPP and, consequently, phytoplankton concentrations. However, we did not address these aspects and further studies are needed to analyse their influence on our parameterisation of phytoconvection.
Besides this the analysis of the new parameterisation revealed some other points requiring improvement. Considering that the MLD is the main controlling factor in both the paramerisation of phytoconvection and in the transition between the winter and summer regimes, the use of a higher vertical resolution between 200 m to 600 m depth would improve its accuracy in spatially explicit, especially during the decline of the MLD. Such improve-ment would also be associated with increased computation times. There is also the need of improving the MLD criterion (9) itself to account for haline stratification. This could be achieved by applying a density-based criterion following the algorithm developed by Kara et al. (2000). The comparison of the MLD criterion applied here with those based on a temperature difference of 0.8 K proposed by Kara et al. (2000) showed that these criteria yield even deeper MLDs which are less consistent with the simulated temperature distribution (see Fig. 5). Hence, further testing is needed to obtain the bestfitting criterion and to include the effect of salinity.
The introduction of the weighting function (11) allowed the smooth transition between winter and summer regimes with deep and shallow mixed layers, respectively. In addition, the weighting function produced improved results in the deep open ocean and shallow shelf regions. However, with respect to the transition from a deep to a shallow mixed layer, further improvement of the weighting between the conventional parameterisation and the parameterisation of phytoconvection is required.
The simulations with the convection model showed that during phases of reduced convective mixing the MLD shows only a slight response, even though these phases last for a couple of days. The simulations also showed that there is a temporal mismatch between the final shutdown of convection and the decline of the MLD. Thus, a possible improvement to the weighting function (11) as well as to the ratio t exp /t orb (see equations (7) and (8)) would be the coupling to the H cml rather than to the MLD.
The vertical exchange coefficients used in some hydrostatic models (as for example the HAMSOM model) in combination with the low vertical resolution are not always sufficient to account for convective mixing of particles. On the other hand, convection-resolving models are not feasible to largescale models due to immense computation times, which represents the major challenge for the further improvement of the here presented parameterisation of phytoconvection. For this purpose, a possible solution could be the coupling of the parameterisation of phytoconvection to the air-sea heat flux, a proxy for convective mixing. Taylor and Ferrari (2011) used the net surface heat flux to show that the shutdown of convection can be a better indicator for the onset of the bloom than the MLD. However, other mechanisms, e.g. eddy-driven stratification (Mahadevan et al., 2012) can lead to stratification without a change in the net surface heat flux, and would, thus, still not be captured by this approach.
Besides these possible and necessary improvements and the further anal-ysis of the underlying assumptions, another important step is the application of the presented parameterisation in a model region where sufficient observation data are available to obtain a significant image of the quality of the parameterisation. Such can be achieved with the help of other largescale, physical-biogeochemical models, used to calculate growth dependency on the local PAR in depth z. The parameterisation of phytoconvection presented here is based on the originally implemented parameterisation of lightdependent primary production, and thus, is expected to be applicable by the vast majority of this type of models.