A large eddy simulation study of the formation of deep chlorophyll/biological maxima in un‐stratified mixed layers: The roles of turbulent mixing and predation pressure

Recent experimental measurements of fluorescence values and turbulent energy dissipation rates, recorded in weakly stratified boundary layers in the open ocean, have highlighted a significant correlation between the formation of deep chlorophyll maxima (DCM) and turbulent mixing. Specifically, the depth of many DCM are observed to lie below, but within about one standard deviation, of the point at which the energy dissipation rate profile reaches its maximum. This correlation of DCM and turbulent mixing is both exciting and curious, as conventional thinking tends to see the latter as a destructive rather than a constructive agent in regards to the formation of deep biological maxima (DBM), for which DCM data is usually interpreted as a proxy. In order to investigate this phenomenon, a three‐dimensional large eddy simulation (LES) of the ocean boundary layer was combined with a generic nutrient‐phytoplankton‐zooplankton (NPZ) type biological model, in order establish what mechanisms might be driving the experimental observations. Simulations of the LES‐NPZ model, based upon various sets of generic biological parameters, demonstrate DCM/DBM formation occurs at normalized depths close to those seen in the experimental observations. The simulations support the hypothesis that the DBM are generated by a combination of zooplankton predation pressure curtailing phytoplankton growth near the surface, and a decline in the strength of the vertical mixing processes advecting nutrient through the boundary layer. In tandem, these produce a region of the water column in which predation pressure is relatively low and nutrient aggregation relatively high, suitable conditions for DBM formation.

The presence of deep (or sub-surface) chlorophyll/biological maxima (DCM/DBM) observed in vertical fluorescence profiles is one of the most ubiquitous features of the world's oceans (Mac ıas et al. 2008), and much research effort has been devoted to understanding the mechanisms behind their formation and dynamics (see Cullen 2015 for a comprehensive review of the subject). Observations of DCM are not just confined to the ocean boundary layer (e.g., Cullen 1982;Estrada et al. 1993;Letelier et al. 2004;Mac ıas et al. 2013), but are a pervasive feature of the limnology of lakes too (Hamilton et al. 2010;White and Matsumoto 2012;Simmonds et al. 2015). Here, the terminology DCM specifies a broad (10-20 m) region of relatively high chlorophyll concentrations (but weak concentration gradients), usually to be found at or somewhat below the mixed layer depth. This should not be confused with the concept of a rapidly varying biological "thin layer," confined to vertical scales of just a few meters (e.g., Dekshenieks et al. 2001;McManus et al. 2003McManus et al. , 2005Benoit-Bird et al. 2009;Durham et al. 2009;Johnston et al. 2009), which are usually found residing within the top 10 m or so of the water column.
A number of different postulates have been advanced to attempt to explain both the formation of DCM, and the reasons as to why they are observed so frequently under many different hydrodynamic conditions. These include the presence of a nutricline, giving rise to increases in phytoplankton growth rates in certain preferential layers of the water column (Simmonds et al. 2015), or the physical accumulation of phytoplankton cells at boundary layer interfaces (Ruiz et al. 2004;Huisman et al. 2006), such as the thermo/ pycnocline (Dekshenieks et al. 2001). One possible mechanism by which aggregations could form at such boundaries would be through the action of intense levels of shear instabilities disrupting the swimming motions and effectively trapping gyrotatic micro-organisms (Durham et al. 2009;Hoecker-Martinez and Smyth 2012). A further alternative, highlighted by the observational and numerical work of Fennel and Boss (2003), emphasizes what they term the "general compensation depth," a level at which the local phytoplankton growth rate is balanced out by losses due to zooplankton grazing and divergences of settling velocity. This emphasis on predation pressure through zooplankton grazing is important because the upper ocean mixing layer most conducive to high levels of phytoplankton growth (it usually encompasses the euphotic zone in which light levels are sufficient for photosynthesis) is turbulent, and turbulence levels significantly influence planktonic predation (e.g., Rothschild and Osborn 1988;MacKenzie and Kiørboe 1995;Lewis 2003;Galbraith et al. 2004;Lewis and Bala 2006 and references therein). Awareness of the role played by turbulence in the formation and sustenance of biological layers has recently been heightened by the publication of new observational evidence. This is derived from data recorded by sophisticated, next generation, instrumentation, which shows that turbulent mixing rates and biological layers are frequently correlated (Wang and Goodman 2010;Mac ıas et al. 2013), highlighting that background turbulence has a potentially creative, as well as a destructive capacity in regard to biological aggregations. The observations of Mac ıas et al. (2013) are particularly striking in this regard, and are discussed in a little more detail in "A brief. . .Mac ıas et al." section. They provide both motivation and empirical support for the subsequent theoretical/computational analysis.
The question that now arises is how to go about investigating the creative and destructive capacity of boundary layer turbulence in regards to DCM/DBM production? The recent increase in instrumentation sophistication has been achieved hand-in-hand with ever increasing advances in both computational speed and storage capacity. This in turn enables one to contemplate the use of much more ambitious and extensive modeling methodologies than have been attempted hitherto. The most common means of modeling planktonic population dynamics has been through the utilization of some form of coupled nutrient-phytoplankton-zooplankton (NPZ) system of differential equations (e.g., Franks 1995;Edwards et al. 2000;Franks 2002). Usually such models contain various simplified descriptions of nutrient uptake, phytoplankton growth, and zooplankton predation built into their equations. In order to study how these processes are effected by background turbulence, one requires a more sophisticated, coupled, biophysical model. Previously, this has been attempted by adding in a vertical eddy diffusivity term into the NPZ system, which in effect creates a one-dimensional (vertical) biological boundary layer model (Baird and Emsley 1999;Denman and Pe na 1999;Flierl and McGillicuddy 2002;Fennel and Boss 2003). However, the use of vertical eddy diffusivities in this way provides only a crude measure of the effects of water motion on the biology, particularly given that turbulence is, by definition, characterized by high levels of vorticity, which is fundamentally three-dimensional in nature (Tennekes and Lumley 1972;McComb 1991;Lesieur 1997). Given the unquestioned advances in computational speeds in recent years, it seems reasonable to try and formulate something more powerful and sophisticated.
In Lewis (2005), the author attempted to devise such a model, building on ideas first formulated in the review article of Denman and Gargett (1995). This was achieved by employing a series of large-eddy simulations (LES) of ocean mixing layers, coupled together with a specially adapted NPZ model (specifically based upon the ideas of Baird and Emsley 1999) to describe the biological evolution in such an environment. LES captures the large scale features of a boundary layer flow extremely well and these are utilized to advect the biological fields in bulk. Although random, on small (planktonic) scales turbulent flows possess certain universal generic features, usually characterized by the local energy dissipation rate e z ð Þ at a specified depth z. These generic features mean that the influence of the physical forcing on the population dynamics will conform to certain statistical norms, which can be estimated. Hence a knowledge of e allows one to parameterize the coupling of the LES to the NPZ, to create a fully integrated, 3D, bio-physical model across all scales. The main details of the adaptations necessary to produce a combined LES-NPZ model of the mixing layer were presented in Lewis (2005). In that work, the behavior of the (initially uniform) biological fields when subjected to a fixed level of wind forcing was investigated. These results were generally satisfactory, but computational limitations restricted the run times to no more than a day or so, severely restricting the interval during which significant heterogeneities in the biological fields might develop. Subsequent advances in computational capacity and processor speed now allow one to carry out runs over periods of days and weeks, giving the opportunity to investigate the potential for the development of aggregations over a number of biological cycles. In conjunction with the new experimental results starting to become available, the LES-NPZ model now provides a means of investigating the drivers of DCM/DBM production and dissipation to a much greater degree of precision than hitherto. The results and analysis of a number of such investigations (presented in "Studies of DCM/DBM formation and characteristics using the LES-NPZ model" section) comprise the main body of this article.
A brief summary of the observations of DCM recorded in the subsurface upper ocean as discussed by Mac ıas et al.
Briefly, the paper of Mac ıas et al. (2013) (denoted by M13 subsequently) describes an extensive series of measurements recorded using a TurboMAP-L fast sampling (512 Hz) probe (Doubell et al. 2009), which is capable of measuring conductivity, sea temperature, vertical shear, and fluorescence on a resolution scale of the order of centimeters. Measurements where m is the kinematic viscosity (which is slightly temperature dependent but a value of m510 26 m 2 s 21 will be assumed throughout this work). M13 do point out that the sampling within the highly turbulent near surface region (where wind stresses and wave effects are at their maximum) may be somewhat problematic, since the TurboMAP-L probe needs to be sinking at an almost constant speed through the water column. Such conditions are only reached at depths 10-15 m below the surface. The prevailing vertical stratification conditions were summarized by means of an average buoyancy frequency N z ð Þ, determined by where q is the water density. Measurements of both fluorescence F z ð Þ (measured in relative fluorescent units) and energy dissipation rate were fitted to standard Gaussian curves to create model depth profiles. Four such profiles, one from each of the locations mentioned above, superimposed on the respective raw data are shown in Fig. 1 (These profiles are reproduced from those of Fig. S1 of M13).
Among the many interesting features of these datasets, the most pertinent to this particular study is the relative position of the deep chlorophyll maximum DCM depth , to the position peak depth at which the energy dissipation rate reaches its recorded maximum. Some 73% of the analyzed e z ð Þ profiles were termed "positive," in the sense that they exhibited a subsurface peak depth . Of these "positive" profiles, two thirds (31 out of 46 profiles) exhibited a DCM depth situated below, but within one standard deviation of peak depth . The average DCM depth for these "matching" (M13) profiles was located some 18 m deeper than the corresponding peak depth values. In other words, the maximum fluorescence measurements were significantly correlated with negative (decaying) e z ð Þ gradients within the upper ocean mixed layer. In the remaining 15 profiles, this correlation was absent. The main difference between these latter mixing layers was their relatively high degree of stratification (N 2 53:21 3 10 23 61:82 3 10 23 s 22 ), compared to the much lower values (N 2 51:1310 23 61:2310 23 s 22 ) associated with the "matching" profiles. These results are important because many theoretical explanations of the formation of heterogeneous phytoplankton concentrations begin with the premise that significant stratification is already present, and such conditions are a prerequisite for biological aggregations (especially thin layer formation, e.g., Durham et al. 2009;Hoecker-Martinez and Smyth 2012). But here is experimental evidence of intense DCM (typical values of the magnitude of DCM max are between two and six times the background) forming in the absence of strong stratification. Usually the presence of a DCM correlates exactly with a corresponding deep biological maximum DBM, although there is some experimental evidence (Longhurst and Harrison 1989;P erez et al. 2006;Cullen 2015) that the latter can be displaced slightly below the former as the chlorophyll-to-carbon ratio increases with depth. However, this displacement usually occurs in relatively stable ocean boundary layers found in tropical waters, conditions which do not apply to any of the locations discussed in M13. Hence for these sites, it is fair to assume that the measured DCM is equivalent to the DBM. So the question arises, what causes the formation of these aggregations in the absence of significant stratification in the upper boundary layer?
M13 postulate that the existence of a DCM may be explained by the juxtaposition of two vertically opposing gradients of resources: light availability from the surface necessary for photosynthesis and vertical mixing of nutrient rich deep waters into the euphotic zone, creating an "optimal window" for the formation of the DCM. This seems a highly plausible hypothesis, albeit one that is likely to be modulated by other factors. M13 highlight the possible role of cells sinking faster out of high turbulence zones and accumulating in low vorticity regions situated below peak depth , a hypothesis supported by the experiments of Ruiz et al. (2004). Another potential regulatory factor is that of planktonic predation (summarized, somewhat loosely, as that of a generic zooplankton species feeding on a generic species of phytoplankton), which is known to be strongly dependent upon the level of background turbulence (Rothschild and Osborn 1988). This is another vertical gradient that augments those governing the availability of light and nutrients. Modern high speed computing resources means it is now quite possible to investigate and test such postulates much more systematically (and cheaply) than ever could be done using field trial data alone (vital though that remains). A necessary pre-requisite is the existence of suitable biophysical models of the upper ocean boundary layer, which incorporate both the effects of wind-driven and surface wave generated turbulence, coupled with the biological drivers of plankton population dynamics. These are just the kind of bio-physical characteristics the LES-NPZ model is able to replicate.

A description of the LES-NPZ model
Conceptually, the biological part of the LES-NPZ model is little different from other three state NPZ models in the DCM are to be found just below z5peak depth , where he z ð Þi310 27 m 2 s 23 reaches its maximum value e max . The measurements are taken from four different sites. Profile 6 is from the Alboran Sea; Profile 21 is from the Strait of Gibraltar; Profile 28 is near the Antarctic Peninsula; Profile 31 is from the North Atlantic. Full details of these sites, along with many further similar profiles are presented in Mac ıas et al. (2013). (Continued on next page).
literature (e.g., Keifer and Atkinson 1984;Fasham et al. 1990;Edwards and Brindley 1996, and previous references cited). Three, nondimensional, scalar fields denoted by N x; t ð Þ5N Ã =N 0 , P x; t ð Þ5P Ã =P 0 and Z x; t ð Þ5Z Ã =Z 0 , representative of nutrient (specifically nitrate), phytoplankton and zooplankton (where N 0 in kg m 23 , P 0 and Z 0 in cells m 23 are suitable reference scales) are assumed to satisfy three advection diffusion equations of the general form

Lewis et al. Formation of deep chlorophyll/biological maxima
Here, D=Dt @=@t1u:r, and u x; t ð Þ5 u; v; w ½ 5 u 1 ; u 2 ; u 3 ½ is the (resolved LES) turbulent velocity field. This provides the first component of the physical coupling to the biology and is derived from a spatially and temporally (over one wave period) averaged version of the full 3D Navier-Stokes equations. First derived by Craik and Leibovich (1976), these consist of the equations of continuity r:u50, momentum and energy Dh Dt 1U S :rh5SGS: In Eqs. 6, 7, f is the Coriolis frequency, U S the Stokes drift velocity, x5r3u the vorticity, g is the acceleration due to gravity, h x; t ð Þ5h r 1h 0 x; t ð Þ the temperature field, q x ð Þ5q 0 1q 0 x ð Þ the fluid density and p S 5p1q 0 2u:U S 1jU S j 2 h i =2 a generalized pressure term. The density is assumed to be proportional to the temperature, such that q 0 =q 0 5h 0 =h r and details of the values of these parameters, including the reference density q 0 and temperature h r , are given in Table 1. The Stokes drift velocity was estimated by assuming that the ocean consists of steady, monochromatic deep-water waves (Philips 1977), which (for convenience) are directed along the x-axis (easterly direction). In which case U S 5 U S e 2kz ; 0; 0 À Á , where U S 5rka 2 , a being the wave amplitude, k the wave-number and r5 ffiffiffiffiffi gk p the wave frequency (see Table 1).
The sub-grid scales (SGSs) used to close the equations are based on a standard Smagorinsky (1963) where m T is an eddy viscosity and Pr is a turbulent Prandtl number. Details of the eddy viscosity formulation in terms of the resolved strain rate S ij and a resolution scale L 0 (set to be 1m throughout) are discussed in Lewis (2005). Note the distinction between q 0 , the unresolved part of a scalar quantity q total 5q1q 0 , and q 00 used later to denote a fluctuation in the LES resolved part q x; t ð Þ5q x; t ð Þ5hqi1q 00 x; t ð Þ, derived from Eqs. 6, 7 (In what follows hqi denotes Þdxdy, the instantaneous horizontal mean of a resolved quantity. Similarly hqi T will denote q z ð Þ5 1 T Ð T 0 hqidt, the time averaged horizontal mean, over a specified interval T). For code verification purposes, solutions of the velocity and temperature (pressure) fields were computed from Eqs. 6, 7 over a domain 120 3 120 m 2 horizontally and to a simulation base depth z S 533 50 ð Þ m, utilizing a basic grid of 40 3 40 3 75 114 ð Þ nodes. This implies a regular resolution scale of Dx5Dy53 m and Dz50:45 m (although the vertical resolution is staggered to give greater resolution near the sea surface). In the course of these investigations, many different turbulent boundary layers were generated. Their characteristics will be distinguished by the values of the Stokes drift velocity U S , and the friction velocity U Ã , which determines the strength of the wind forcing applied at the surface m T @u @z Here, s is the surface wind stress, which was varied between simulations over a range of 2:25310 23 -25:0310 23 kg m 21 s 22 . This implies values of U Ã lying between 1:5310 23 m s 21 and 5:0310 23 m s 21 , roughly equivalent to windspeeds U 10 51:2-4:0 m s 21 at a height of 10 m. The corresponding values of U S are given in Table 1. All the boundary layers were made slightly convective, with a turbulent buoyancy flux of w 0 h 0 521:2310 26 K m s 21 applied at the surface. Other boundary conditions imposed on the flow are horizontal periodicity, w50 at z50 and no slip at z52z S . Typically, the various boundary layers were spun up from rest for a period s spin $ 60; 000 s until a quasi-equilibrium state was reached, before any biological fields were added. This marks time zero for the simulations proper. The physical characteristics governing the generation of a particular simulation will be denoted by its U Ã ; U S ð Þnumber. So the notation U Ã ; U S ð Þ 5 3:5; 3:9 ð Þ specifies a simulation in which U Ã 53:5310 23 m s 21 and U S 53:9310 22 m s 21 ; respectively. The corresponding Langmuir number La5 ffiffiffiffiffiffiffiffiffiffiffiffiffi U Ã =U S p , as defined by McWilliams et al. (1997), can also be used to classify the simulations.
The performance of the LES code was extensively tested by comparing its output with that of a similar code employed by McWilliams et al. (1997) to study Langmuir turbulence in the upper ocean boundary layer. Comparison profiles of horizontal currents, shear stresses, velocity variances, and energy dissipation rates (Lewis 2005) showed excellent agreement across the mixing layer z 2 2z ML ; 0 ½ . Here, z ML represents the mixing layer depth (e.g., de Boyer Mont egut et al. 2004), which varies according to both the strength of the wind forcing and the local hydrography (density and/or temperature differences, usually measured from the surface). The main difference between the two codes is the incorporation by McWilliams et al. (1997) of a stably stratified region below z ML extending to a depth of z % 23z ML . By contrast, the LES-NPZ model does not feature any such stratification. This raises a problem with the use of the terminology "mixing-layer depth" when applied to these simulations, since the hydrographical features that help regulate its position are absent. So for this work, an alternative concept of turbulent or Ekman depth, will also be used. This depends purely on the mixing properties of the flow and is defined by Coleman et al. (1990) to be z TD 5U Ã =f . In practice, z TD will be larger than z ML , since the latter is constrained by the hydrographic features of the flow not incorporated into these simulations (Pearson et al. 2015). A rough equivalence would be closer to z ML $ z TD =2. From a review of the relevant literature, Garratt (1992) suggests z ML $ 0:2-0:4z TD for weakly stratified layers, but the inclusion of the Stokes drift term in Eq. 6 means that the mixing layer of these simulations is certain to be somewhat deeper than this. Since the velocity field is virtually zero below z TD , while the majority of M13's measurements of DCM are found within weakly stratified mixed layers at depths significantly above z TD (the turbulent mixing, although declining, is still a prevalent feature in Fig. 1), the lack of a stratified base to the computational domain is not a serious issue for this work.
Figures 2-4 show some typical flow parameters derived from samples of these simulations. Figure 2 shows the depth dependence of the mean (LES) velocities hui T and hvi T derived from a U Ã ; U S ð Þ5 3:5; 3:9 ð Þ simulation after a time T5T sim $ 20-25 d. Typically, except very near the surface, the mean flow is directed in a south-southwest direction, with both hui T and hvi T negative. The effect of the Stokes drift term is to deflect the mean current anticlockwise-away from the easterly wind direction-considerably beyond the 45 predicted by classical Ekman boundary layer theory (Lewis and Belcher 2004). Figure 3A shows a comparison of the mean velocity variance profiles hu 00 2 i T ; hv 00 2 i T ; and hw 2 i T taken from the same simulation. (N.B. ww 00 since hwi50.) A feature of these results is the relatively large scale of the vertical hw 2 i T , as opposed to the horizontal hu 00 2 i T 1hv 00 2 i T turbulent mixing. Physically this is a manifestation of the formation of Langmuir cells, a series of counter rotating vortices aligned with the wind direction, driven by the interaction of the wave field stretching and tilting the vertical vorticity field (Teixeira and Belcher 2002). From a biological point of view, this enhanced mixing helps to ensure any nutrient resources, drawn from the deep ocean by large scale upwelling events, will be distributed quickly and uniformly throughout the mixing layer.
The other key physical parameter which encapsulates the mixing is the mean energy dissipation rate e z ð Þ highlighted earlier (q.v. Eq. 1 and Fig. 1). Figure 4 shows some profiles of e z ð Þ, derived by balancing the (steady) resolved scale turbulent kinetic energy budget (McWilliams et al. 1997;Lewis 2005, Eq. 11;Polten and Belcher 2007), extracted from the 2:0; 2:2 ð Þ , 3:5; 3:9 ð Þ ; and 5:0; 5:6 ð Þsimulation datasets over a period of T % 1 d. This involves combining terms such as Stokes, shear and buoyancy production, pressure working, turbulent transport, and SGS dissipation. Near the surface e z ð Þ shows significant variation, by roughly a factor of five in magnitude as U Ã increases, but at greater depths the differences are far less evident. Notice too that the e z ð Þ profiles often exhibit small secondary maxima, usually lying between about 5 m and 15 m below the surface. Such secondary maxima are a feature of these wind and wave driven Langmuir boundary layers. They can be explained by a comparison with the velocity variance profiles of Fig. 3A. While the horizontal hu 00 2 i T 1hv 00 2 i T mixing component declines uniformly, this can sometimes be offset by vertical hw 2 i T component which initially increases near the surface, giving rise to the small secondary maxima in e z ð Þ. Comparing the profiles derived from these numerical simulations with the experimental measurements of e z ð Þ recorded by M13 (see Fig. 1), one can see that the values all lie in a broadly similar 10 27 -10 26 m 2 s 23 range. The main distinction is that for the numerical simulations the peak dissipation rate lies at the surface, where the wind and wave forcing are at their most intense, while in the experimental results the peak depth lies somewhat deeper. Careful inspection of all the low stratification profiles recorded off Spain and Gibraltar (e.g., Profiles 6 and 21 in Fig. 1) shows that the mean peak depth occurs at a relatively shallow 28.9 m (M13 Table S1). In addition, most of these profiles exhibit dissipation values between 1 3 10 27 m 2 s 23 and 2 3 10 26 m 2 s 23 at depths 10-15 m below the surface, similar to those seen in the numerical simulations.
The profiles recorded at the North Atlantic and Antarctic Peninsula sites (e.g., Profiles 28 and 31 in Fig. 1) are somewhat different, exhibiting deeper values of peak depth $ 70 m. The probable reason for these large sub-surface peaks is an increase in the energy dissipation rate at the mixed layer/thermocline boundary. Observations of increased energy dissipation rates within the thermocline have been recorded (e.g., Moum and Osborn 1986) and are believed to be result of the action of internal wave scattering due to buoyancy fluctuations (Gregg 1989). However, the magnitude of this increased mixing generally lies within the range 10 210 -10 29 m 2 s 23 , which is rather too small to account for peaks shown in Fig. 1 and there may be other factors at work (see "Sensitivity to wind forcing" section). In any event, this should not distract from the main feature of all the profiles, namely the pattern of DCM formation at depths where e z ð Þ is in sharp decline, which is motivation for this work. The biological source/sink terms making up the righthand sides of Eqs. 3-5 are inspired by the work of Baird and Emsley (1999) and discussed in detail in Lewis (2005), so only a brief resum e is included here. Table 2 lists the main numerical values of the various biological parameters that appear in these terms. The nutrient uptake term in Eq. 3 is given by provided Eq. 10 is positive, otherwise it is set to zero. This reflects the balance between nutrient uptake by means of turbulent diffusion and the limitations brought about by the nitrate storage capacity of an average cell. Here, r p is a typical phytoplankton cell radius, Sh, a nondimensional turbulent Sherwood number, D N ; the molecular diffusivity of nitrate, and R N 0 z ð Þ=R max N is the ratio of the nitrate storage capacity of the cell to its maximum potential storage capacity. The value of R N 0 z ð Þ, assuming a background ambient nitrate concentration level N 0 set to be 2:8310 25 kg m 23 (Fasham et al. 1990), can be established using a mass balance equation (Baird et al. 2001), as discussed in Lewis (2005). The R max N parameter was set to a value of 3s N where s N (kg cell 21 ) is a nitrate stoichiometry coefficient quantifying the minimal amount of nitrate needed for a cell to be viable. Whenever the background concentration N x; t ð Þ greatly exceeds N 0 then nutrient satiation sets in, at which point diffusion into the cell ceases. For the simulations results presented in the "Studies. . . model" section, the value of R N 0 z ð Þ=R max N % 0:5, which means that uptake ceases when the background nutrient concentration N x; t ð Þ > 2. The physical forcing influences the scale of Eq. 10 in two ways. Directly, through the dependence of Sh on e, and indirectly through the effects of the turbulent mixing on the nutrient distribution. In practice, the direct influence of e on Sh is weak, since r P is usually less than the Kolmogorov length scale m 3 =e À Á 1=4 , meaning that Sh varies little from unity across the boundary layer. The indirect influence of e on Eq. 10, as manifested through the nutrient distribution is far more important, and will form part of the investigation discussed later on. The nitrate recycled and phytoplankton growth terms employed in Eqs. 3, 4 are essentially analogues of each other and depend significantly on the scale of Eq. 10. In the model, they are given by P growth5b E min 1; Equation 12 encapsulates the fact that a phytoplankton species can potentially reproduce at its maximum growth rate l max P under ideal conditions; but this is regulated by the  available light intensity for photosynthesis, which is assumed to decay exponentially with depth (a is the light attenuation coefficient for water), and the cellular reserves of nitrate currently available (see Baird and Emsley 1999 for more on the details underpinning these terms). Other factors can also inhibit phytoplankton growth efficiency and their effects are represented by the dimensionless parameter b E 2 0; 1 ½ . These growth inefficiencies (cell death and decay) usually lead to a recycling of nitrate back into the water column, a process which is summarized by Eq. 11. The usual initial condition employed in these simulations will be to fix N x; 0 ð Þ51, in which case the recycling term is about ten times smaller than the uptake term (Eq. 10).
The loss of phytoplankton to zooplankton grazing and the zooplankton growth rate terms used in Eqs. 4, 5 are also analogous. In the model, they appear as Here, I R; T R ; e; r P ð Þis representative of a predation rate integral (units m 3 s 21 ) defined explicitly in Lewis and Pedley (2001, Eqs. 16, 17), evaluated over a range 0; R ½ , where R is the contact radius of the zooplankton predator. This parameter is the maximum distance over which it can perceive its prey (R $ 13 10 23 -40310 23 m) and in Lewis and Pedley (2001), it is assumed to be spherically symmetric for mathematical simplicity. In reality, most predators possess a much narrower conical perception field, the ramifications of which are discussed in Lewis (2003), Bala (2006, 2008). Loosely I R; T R ; e; r Z ð Þestimates the predation rate of the predator from its encounter rate (Rothschild and Osborn 1988), times its capture efficiency. The latter is governed by parameters such as a predator's reaction time T R , its swimming or pursuit speed r Z and level of turbulence summarized by e. Turbulence always enhances encounter rates by advecting more prey particles into the predator's vicinity, but in certain instances can lead to the suppression of the predation rate, because it makes the act of actually capturing prey more problematic. In this work, the zooplankton predators are relatively small r Z $ 10 25 210 24 m, and their corresponding swimming speeds r Z $ 5310 24 2 20310 24 m s 21 are an order of magnitude smaller than the mixing velocity scales found in the upper boundary layer (Fig. 2). Consequently, swimming is not included directly into the advection terms in Eq. 5 because it would be too small to influence the spatial distribution of the zooplankton. Potentially, swimming could make a significant impact at depths close to the mixing layer boundary where e falls off rapidly and DBM are observed. This is particularly true of the predation term, which falls to zero in the absence of any advective or swimming motions. On balance, this would a fairly unlikely scenario, but for illustrative purposes nonswimming predators will feature in some of the simulations (cf. Figs. 14, 16) carried out here. The specific ideas as to how these competing factors are formulated into the integrand of I R; T R ; e; r Z ð Þ , are discussed in Lewis and Pedley (2001). Figure 5A shows some profiles of the predation rate I R; T R ; e; r Z ð Þfor a relatively efficient T R 55 s ð Þnonswimming predator, possessing a contact radius R52310 23 m, over a range of windspeeds (U Ã 52:0 310 23 -5:0310 23 m s 21 Þ. Generally, an efficient predator benefits from increasing levels of turbulent mixing, since its encounter rate goes up leading to more prey captures. This is reflected in Fig. 5A which exhibits a steady increase of predation rate with windspeed. The benefit is particularly marked near the surface where e is at its largest, q.v. Fig. 4. By contrast, increased turbulent mixing can be detrimental to the predation rate of an inefficient predator T R 515 s ð Þ , such as the one illustrated by the profiles in Fig.  5B. In these instances, the comparatively high relative velocities found near the surface means the predator has to react faster than T R in order to capture prey moving into its vicinity. This it is unable to do, leading to a drop in its overall predation rate. The predation rate determines both the decline in the general phytoplankton population through grazing (Eq. 13), and the zooplankton growth rate (Eq. 14). The latter is regulated by the yield Y of new predator cells per prey cell captured, and is also restricted never to exceed a fixed maximum zooplankton growth rate l max Z , theoretically attainable under ideal conditions.
The biological part of the model is closed by assuming that the zooplankton growth rate is limited by a simple constant mortality term of the form The mortality rate l Zdeath is a purely biological parameter independent of e. Its importance stems from the fact that it regulates the period of oscillations of the planktonic P=Z populations within the model (see "Model . . . dynamics" section). The diffusion coefficients in Eqs. 3-5 are calculated from the eddy viscosity by means of a turbulent Schmidt number The Schmidt number cannot be too large as it is regulated by the LES resolution scale and one cannot expect to resolve the scalar fields down to a finer level than that can be achieved for the velocity fields. A proposal made by Sullivan et al. (1994) and adopted in Lewis (2005) was to prescribe Sc turbulent in the form where D 3 5 3Dx=2 ð Þ3Dy=2 ð Þ Dz. For the simulations discussed here Sc turbulent % 1=2, a value compatible with L 0 scale which prescribes the resolution of the LES velocity field.
Equations 3-5 are solved subject to certain prescribed boundary conditions on the biological fields. Horizontally, periodic boundary conditions were imposed. Vertically, zero surface flux conditions are (usually) imposed for each field, while at z52z S certain prescribed fluxes into the simulation domain were enforced. For the majority of the simulations presented here, the vertical boundary conditions satisfied @C @z z50 50; u T Sc turbulent @C @z z52zs 5hwCi; C5N; P; Z ð Þ : (18) In the case of the planktonic fields, C5P and Z, deceased cells were assumed to sink, under gravity, out of the simulation domain at a rate governed by a settling velocity U sink (Lewis 2005). This fixes hwCi5U sink / dead C C 2z S , where / dead C is the (small) proportion of deceased cells out of the total population. The boundary conditions imposed on the nitrate field were more flexible. In general wind, tidal and other driving forces instigate upwelling motions that circulate nutrient rich water from the deep ocean through the mixing layer. So one would expect nitrate levels in the mixing layer to remain relatively constant, as losses through uptake by the phytoplankton for photosynthesis will be quickly replenished via this source. Hence a positive flux of nitrate into the simulation domain was invariably imposed. Williams and Follows (1998) suggest this background flux should be $ 2 3 10 28 mol N m 22 s 21 , which is roughly equivalent to hwN Ã i522:8 3 10 210 kg m 22 s 21 . This figure was taken as representative of the background nutrient replenishment flux from the deep ocean into the boundary layer. However, in instances where z S ) z TD , introducing replenishment nutrient at the base is not very useful, since there is little or no mixing to advect it through the boundary layer. To overcome this problem, simulations were also conducted in which the replenishment flux was imposed at the surface (z50) instead. Such a surface nutrient surge could be driven by a significant river run off event, iron seeding experiments or even ash from a volcanic eruption (Frogner et al. 2001). This allows one to investigate what changes (if any) occur to the DCM/DBM characteristics when the replenishment nutrient is added at alternative depths subject to radically different turbulent mixing regimes.

Model verification and the influence of boundary layer structure on the biological dynamics
The full version of the LES-NPZ model seeks to summarize the effects of a number of different physical and biological drivers on the planktonic populations. It depends upon many different parameters, not all of which are critical in determining the evolution of the biology. To determine which terms are important and explore a suitable parameters space for the biological variables, a greatly simplified version of the NPZ model described above was developed with all the advection/ diffusion terms and boundary layer structure removed. This simplified NPZ model reduces to the formulism where H x f g is the Heaviside step function (H x f g51 if x ! 0 and zero if x < 0) and the other constants are given by These terms no longer exhibit any dependency on depth z or turbulent mixing e, which are features of the full LES-NPZ model. In this reduced form, Eqs. 19, 20 are sufficiently simplified to allow some qualitative analysis of their behavior. Assuming hN t ð Þ < 1 and cP t ð Þ < l max Z , (the most common scenario) the equations possess a single co-existence equilibrium point in positive phase space, situated at Employing some fairly typical biological parameters, as listed in Table 2, one finds h 21 > 10ed 21 , hence the approximation made in Eq. 21. Linearized stability analysis of the community matrix of partial derivatives at this equilibrium point yields three eigenvalues of the form These eigenvalues give rise to three solution branches near N EQ ; P EQ ; Z EQ À Á , one of which is an exponentially increasing on a time scale of s exp 5c=dl Zdeath , while the remaining two associated branches give rise to oscillatory solutions for the P; Z ð Þ fields on a timescale s oscil 52p= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi al Zdeath p . Based upon the biological parameters specified in Table 2 (specifically for those used in Fig. 6), one finds that s oscil % 10:5 d while s exp % 43 d. So this analysis suggests that if the initial populations happen to lie close to N EQ ; P EQ ; Z EQ À Á % 1:56; 0:43; 0:85 ð Þ , then P; Z ð Þ will oscillate regularly about their equilibrium values at a relatively rapid frequency $1=s oscil , while growing relatively slowly over the longer time scale s exp . However, this conclusion is slightly deceptive because the nutrient uptake term will shut down if ever N gets too large and cells become satiated (when hN > 1 in Eq. 19a). In such circumstances, the predicted growth branch actually manifests itself as a decay branch, over the same long timescale s exp . This means that the P; Z ð Þ populations simply oscillate regularly, while the nutrient concentration exhibits a net depletion, brought about by funding of the new P growth. Qualitatively, this kind of behavior is illustrated in Fig. 6, which shows the evolution, over a period of 25 d, of the normalized Phytoplankton-Zooplankton-Nitrate fields, starting from initial values of N; P; Z ð Þ t50 5 1:0; 0:5; 0:5 ð Þ . Oscillations in the P; Z ð Þ fields over a timescale close to the predicted s oscil % 10:5 d are apparent, while the N field decays relatively slowly from its initial value, since there is no mechanism for nutrient replenishment in the simplified NPZ model. Note too how the peaks in the Z field always lag a few days behind those of the phytoplankton concentration. This kind of behavior is generic to other NPZ models, e.g., Baird and Emsley (1999), Edwards et al. (2000), and Franks (2002).
In the simulations of the full LES-NPZ model which follow, the evolution of the P and Z-fields exhibits qualitatively similar oscillatory behavior to that illustrated in Fig. 6. The period of the oscillations in the mixing layer always lies close to the value set by s oscil , except in those instances when the initial conditions P; Z ð Þ t50 differ significantly from the equilibrium values P EQ ; Z EQ À Á . The default values of P; Z ð Þ t50 50:5 usually employed are close enough to equilibrium for oscillatory behavior to commence immediately. The incorporation of turbulent mixing brings about some quantitative changes in the biological dynamics. First, the peak concentrations are somewhat reduced and smoothed out compared to those shown in Fig. 6 (typically by around 40% or so). The biological parameters in the simplified model are set appropriately for conditions at the surface z50, but in reality no phytoplanktonic cell will reside so close to the surface all the time and as a result its growth will be retarded by a lack of light when carried deeper into the boundary layer. This is counter-acted, to some extent, by the reduction in zooplankton predation pressure as one descends deeper into the boundary layer. The combination of these effects reduces the peak concentration scales. Second, as will be seen in many of the simulations, the turbulent mixing does not extend uniformly across the entire simulation domain. This means a cell initially situated somewhere near the top of the mixing layer has a higher probability of residing  Table 2. near the surface for longer periods and replicating more frequently than a cell initially residing below the mixing depth (Since it experiences greater exposure to elevated light levels). This effect gradually manifests itself in a visible lag in the intervals separating successive concentration peaks at different depths. The nutrient flux boundary condition (Eq. 18) is sufficient to replenish nutrient losses due to phytoplankton growth, so the decline in nutrient concentration exhibited in Fig. 6 is not a feature of the full LES-NPZ model.

Studies of DCM/DBM formation and characteristics using the LES-NPZ model
Before looking in more detail at those bio/physical parameters which most influence the formation of DCM/DBM, it is worth pointing out some generic features common to all the LES-NPZ simulations results reported here. Each boundary layer is generated by means of the fixed wind forcing boundary condition (Eq. 9), summarized by its U Ã ; U S ð Þparameters. Typically a boundary layer is spun up from rest and the respective U Ã ; U S ð Þ forcing applied so that it relaxes into a state of quasi equilibrium, a process which takes about s spin $ 60; 000 s. It is at this point the mean e z ð Þ profiles are estimated and these appear alongside the biological profiles in the figures that follow. Knowledge of e z ð Þ allows one to calculate the biophysical coupling terms that appear in Eqs. 3-5. The whole process is then repeated, but this time the onset of quasi-equilibrium is the signal for the biological fields to be introduced into the boundary layer. Biological evolution (usually) commences from an initial distribution N; P; Z ð Þ t50 5 1:0; 0:5; 0:5 ð Þ , uniformly applied both horizontally and vertically. The biological fields are then allowed to evolve for a period T sim $ 20-25 d, in accordance with the biological parameterizations set out in Table 2.
The choice of simulation time period T sim $ 20-25 d is driven primarily by biological considerations. One needs to allow the biological fields to replicate through at least a couple of reproductive cycles. These values of T sim are somewhat larger than s pred the predictability timescale of the LES. The latter is the time by which imperfections brought about by not resolving down to the smallest scales propagate through the eddy hierarchy to produce significant contamination of the large scale motions (see Lesieur 1997, Chapter XI, XII). For this work, s pred % 25T E 525 3 large eddy turnover scale % 1:5 d. This does not invalidate the methodology, because a LES computed over times T sim > s pred is still a statistically realistic representation of an actual flow, but one that is in the process of being advected as a complete body of water from its original position (as would happen under the action of a large scale current encompassing a much greater volume than the typical 120 3 120 3 33 or 50 m 3 domains used here). It is important that the flow statistics generated by the LES remain roughly stable over a simulation time period. Figure 3B shows the evolution of the mean velocity variance hu 00 2 i T profile over a T sim 523 d period, taken from a representative U Ã ; U S ð Þ5 5:0; 5:6 ð Þ simulation. As one can see the variance remains almost constant throughout, and this scenario applies for both the hv 00 2 i T and hw 2 i T variances too. So the basic statistics are sufficiently robust for meaningful conclusions to be drawn.

Sensitivity to wind forcing
The first set of simulations were designed to investigate the sensitivity of DCM/DBM formation to wind forcing. To achieve this, three different wind driven boundary layers (each with a simulation depth z S 5250 m) were spun up, with parameter settings U Ã ; U S ð Þ5 2:0; 2:2 ð Þ ; 3:0; 3:3 ð Þ , and 4:0; 4:4 ð Þ . These values are representative of low, intermediate, and strong wind forcing regimes, respectively (the terminology "strong" is only comparative, since it represents a windspeed of only around 4 m s 21 ). Figure 7A-C show the energy dissipation rate profiles for each of the three boundary layers in turn, alongside the corresponding profiles of zooplankton predation. Derived from the initial spin up run, these average profiles remain fixed throughout the subsequent biological simulation. The dissipation profiles show a progressive increase in the near surface maximum dissipation rates (0.65 310 27 -2.0310 27 m 2 s 23 ) with wind forcing, and a corresponding deepening of the mixing layer. In each case, the dissipation rate falls to (almost) zero around about half the turbulent depth % z TD =2, equivalent to 10 m, 15 m, and 20 m, respectively. The predation rate profiles I R; T R ; e; r Z ð Þ are calculated for a relatively large predator length scale r Z 53310 24 m À ), possessing a spherical perception field of radius R51:3310 22 m, with an average swimming speed r Z 52310 24 m s 21 , reacting in a time T R 55 s (see Table 2). These profiles are all fairly similar near the surface (as T R is only moderately fast, the predator cannot benefit from the increase in prey contacts brought about by the highest surface dissipation rates), falling away at different rates, before levelling off $ 10 27 m 3 s 21 near z TD =2. Since e z < z TD =2 ð Þ%0 and the flow is relatively quiescent, the predator is forced to rely on its swimming capabilities alone in order to find prey. This means the predation rate is almost constant for z < z TD =2. Figure 8A-C show the corresponding evolution of the biological N; P; Z ð Þprofiles over the simulation time T sim , assuming a maximum (surface) phytoplankton growth rate of l max P 55310 25 s 21 . In these simulations, nutrient losses through biological growth were compensated for by a nutrient flux hwN Ã i into the boundary layer through the base z S of the simulation domain. Since losses from the former are easily outweighed by gains from the latter, the effect over time is to create a nutricline, with significantly higher nutrient concentrations at the base of the layer compared to the surface. Consider the U Ã ; U S ð Þ5 2:0; 2:2 ð Þcase Fig. 8A first. The mixing layer of z ML % z TD =2510 m is relatively shallow, insufficient for any of the extra nutrient added at the base to  Þforcing. Biological data as given in Table 2, incorporating the large zooplankton predator.
be vertically mixed into the surface region. So it simply accumulates, uselessly, in the lower layers. However, this is of little import since the initial nutrient level N51:0 is already sufficient to support immediate P growth. However, this growth is not uniformly distributed across the boundary layer. Instead a DBM forms at a depth z58-15 m, at or just below the point z ML % z TD =2 where e z ð Þ % 0. The position and form of the simulated DBM is strikingly similar to those recorded in the experimental datasets shown in Fig. 1 (which are for DCM). This illustrates the importance of M13's results, since most other observational reports on DCM do not include any corresponding measure of turbulent mixing. This numerical modeling work strongly supports M13's observations that the two are indeed correlated. Notice too that DBM is not brought about by the presence of the nutricline, because that has not had time to form as yet. The numerical DBM is a transient feature, lasting about 5 d or so because, just as in the simplified model illustrated in Fig. 6, the growth in the P concentration promotes a corresponding surge in the zooplankton concentration. A numerical DZM forms at a similar depth but somewhat after the DBM, consuming the latter before itself dies off. This restores the biological fields to a uniform (not the initial) state after about 20 d. One would anticipate the cycle would be repeated roughly in accordance with the dynamical properties discussed for the simplified model (for these parameters s oscil % 12 d at a depth of z % 7 m), although this has not happened by the end of the simulation. For the U Ã ; U S ð Þ5 3:0; 3:3 ð Þsimulation shown in Fig. 8B, the results are broadly similar. Again the z ML % z TD =2515 m depth is still too shallow for much of the replenishment nutrient added at z S 550 m to be mixed into the surface region. But with initial nutrient levels already high across the simulation domain, a DBM still forms relatively quickly. As before the DBM forms just below z TD =2, between z517-23 m, persists for about 5 d, before being consumed by the numerical DZM, which itself dies out through lack of food after about 15 d. Notice the initial DBM is thicker than in Fig. 8A but less intense. This fits with statistical analysis of M13 (specifically Fig. 4 of that paper) that within the DCM depth 60:5DCM thick range, the highest fluorescence values are associated with smaller values of e z ð Þ. However, for the strongest wind forcing case (Fig.  8C) the results are quite different. In this instance no DBM forms. Instead the P profile grows almost uniformly across the whole of the simulation domain, before being consumed by the growth in zooplankton it has stimulated.
It is worth looking more closely at why no DBM forms for the 4:0; 4:4 ð Þcase, when based on the other less windy simulations, one might expect a somewhat less intense maxima to occur somewhere about z525 m. There are two possible reasons. First, the e z ð Þ profiles derived turbulent kinetic energy budget equation, have a tendency to underestimate the extent of the vertical mixing lower down in the boundary layer. This can be seen by examining the horizontally averaged vertically velocity variance profiles hw 2 i T for the three simulation regimes (Fig. 9A-C). For the low wind case 2:0; 2:2 ð Þ , the hw 2 i T profile closely matches the corresponding e z ð Þ profile (cf. Figs. 7A, 9A), but for the higher wind cases the hw 2 i T profile decays more slowly and penetrates to a deeper level than the corresponding e z ð Þ profiles. This behavior is also a feature of other LES ocean boundary layer codes (e.g., McWilliams et al. 1997;Pearson et al. 2015) which incorporate Langmuir circulations. It is also evident from M13 experimental measurements of e z ð Þ (Fig. 1) that energy dissipation does not abruptly cease at one particular depth, but rather there is evidence of intermittent turbulent bursts extending below peak depth 2peak thick . Polten and Belcher (2007) point out that in contrast to classic shear driven boundary layers, inclusion of significant wave effects characterized by the Stokes drift term U S , results in increased vertical transport associated with the formation of downwelling jets carrying fluid down as far as the turbulent depth z TD . This feature is probably a factor underlying the strong sub-surface peaks observed in the e z ð Þ profiles recorded at the North Atlantic/Antarctic sites. For the strong wind/wave 4:0; 4:4 ð Þsimulation, this enhanced vertical mixing penetrates to z TD 540 m, preventing the formation of a DBM in this instance. The other possible reason is that in this instance the setting of z S 550 m is too small, and if the simulation domain were to be extended to say z S 575 m, it would be large enough facilitate the generation of a weak DBM. M13's observations show that DBM can occur at depths down to 100 m or so, in boundary layers subject to genuinely strong U Ã > 5:0310 23 m s 21 wind forcing conditions. Figure 10A-C shows the evolution of the biological N; P; Z ð Þ profiles for the same three boundary layers as before, only this time the replenishing nutrient flux is applied at the surface rather than the base. So unlike in the previous examples, this extra nutrient is mixed throughout the mixing layer, eventually reaching depths close to z TD for the 4:0; 4:4 ð Þsimulation. Although a surface flux is a somewhat unrepresentative of how nutrient replenishment typically occurs within ocean mixing layers, it provides a means, within the model constraints, of making the extra nourishment readily available to the biological populations in those regions of the water column where DBM formed previously. For the low and intermediate wind simulations 2:0; 2:2 ð Þand 3:0; 3:3 ð Þ , the initial behavior is much as before, with transient DBM forming at around z ML % z TD =2 (the extra surface nutrient means the DBM form slightly closer to the surface and are more intense). However, the main effect of adding replenishment nutrient at the surface is to stimulate the appearance of a secondary DBM, about 15 d after the primary DBM was consumed by the zooplankton population. The secondary DBM develop at the same depths as the primary and last about the same time, but are slightly less intense. For the strong wind 4:0; 4:4 ð Þ simulation, the extra vertical mixing prevents the formation of either a primary or secondary DBM, much as before. Two points come out of this. First, Þforcing. Biological data as given in Table 2, incorporating the large zooplankton predator.
DBM formation is clearly a fairly robust process, which operates independently of changes in the nutrient distribution within the boundary layer. Provided enough nutrient is present initially and the mixing is not too strong, DBM are very likely to form. Second, if background nutrient levels remain fairly high through some effective means of replenishment as in Fig. 10, there is a tendency for the characteristic biological timescales to manifest themselves over and above the homogenizing tendencies of the physical boundary layer drivers. Hence the formation of the secondary DBM, seen in the low wind simulations, occurs after about 15 d, close to the predicted s oscil % 12 d timescale derived from the biological parameters alone. By contrast, secondary DBM had not formed by the end of the base flux simulations, because not enough nutrient was present to stimulate them.

Sensitivity to predation pressure and nutrient initial conditions
The simulations presented in the previous section assumed a relatively large (r Z 53310 24 mÞ zooplankton predator, with a correspondingly large contact radius R51:23 10 22 m and swimming capabilities r Z 52310 24 m s 21 . It is interesting to investigate what happens if the predator were made somewhat smaller, reduced in size to r Z 55310 25 m, with the other parameters rescaled accordingly (see Table 2, e.g., r Z 55310 25 m s 21 and R52310 23 m). The effect of these basic size reductions is to dramatically reduce the predation rate for an individual predator. This is illustrated in Fig. 11, which shows the corresponding predation rate calculated from a e z ð Þ profile taken from a U Ã ; U S ð Þ5 3:0; 3:3 ð Þsimulation (with z S 5233 m rather than z S 5250 m). Compared Þforcing. An inward nutrient flux boundary condition is implemented at the base z s 5233 m of the simulation domain. Biological data as given in Table 2, incorporating the small zooplankton predator.
to the corresponding profile for the large predator (Fig. 7B), the individual capture cross section has fallen by a factor of $ 200. However, smaller predators tend to be more numerous than larger ones. In these illustrative examples, the background concentration for the large predator was set to 10 2 cells m 23 , increasing to 2310 4 cells m 23 for the smaller predator. This effectively offsets the drop in capture cross section, meaning the overall predation pressure on the P field is little changed. Figure 12 shows the evolution of the biological N; P; Z ð Þ profiles based upon the small zooplankton predator with replenishment nutrient added at the base. Compared to the corresponding evolution profiles U Ã ; U S ð Þ5 3:0; 3:3 ð Þ for the large predator, shown in Fig. 8B, a couple of distinctive features are obvious. On this occasion, the initial DBM forms at just below z5z TD =2 between about 15 m and 20 m much as before, but then two secondary DBM evolve, after about 10 d and 20 d, respectively. The reason for this much faster evolution lies in the fact that the death rate of the small predator has increased by a factor of 10=3 compared to that of the large predator (smaller organisms tend to have shorter life-spans). From the analysis of "Model . . . dynamics" section, this reduces s oscil 52p= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi al Zdeath p to about 6:5 d, leading to the production of two secondary DBM over the T sim interval. The other distinctive feature is that the secondary DBM are almost as intense as the primary and form a little deeper in the water column. This can be put down to the fact that the simulation domain is somewhat shallower than previously, which means that it is easier for nutrient added at the base to reach the regions where DBM formation occurs. So it is more easily accessible for the phytoplankton, producing a strong response in the secondary growth phases. Looking at the corresponding results for a surface nutrient flux in Fig. 13 (cf. Fig. 10B for the large predator), the results are broadly similar to the base flux results, with the formation of two secondary DBM, although on this occasion they are slightly less intense than the primary. Since nutrient entering at the surface is mixed more thoroughly throughout the mixing layer, the background concentration never reaches the levels seen in the base flux results, reducing the intensity of the secondary growth surge. Notice too that the mixing has the effect of gradually merging the secondary maxima together, so the DBM becomes a more permanent feature of the boundary layer. Overall, these results show that phytoplankton populations residing in the turbulent boundary layer will tend to evolve to form DBM a little below z TD =2, provided the boundary layer contains enough nutrient to stimulate them. These DBM are robust features of the biological dynamics, and form independently of the finer details of the predation pressures imposed by different zooplankton species. So observations of DBM formation in real biophysical boundary layers, which may contain a myriad of co-existing/competing planktonic species, should not be surprising. Þforcing. An inward nutrient flux boundary condition is implemented at the surface z50 m of the boundary layer. Biological data as given in Table 2, incorporating the small zooplankton predator.
These conclusions lead to the question what would happen if the initial nutrient conditions were not conducive to DBM formation? To this end, a simulation was carried out with N; P; Z ð Þ t50 5 0:1; 0:5; 0:5 ð Þ , to see what effect a very small initial nutrient concentration would have on the development of the biological fields. In this simulation, the wind forcing was relatively low U Ã ; U S ð Þ5 2:5; 2:8 ð Þ , and nutrient replenishment was added through the surface. The phytoplankton field was subject to a predation pressure term similar to that shown in Fig. 11, except that in this case the swimming contribution was switched off r Z 50:0 m s 21 . This was compensated for by a reduction in the reaction time to T R 51 s. Since the predator can no longer swim, it means that the phytoplankton field is not actually subject to any predation pressure below the turbulent depth z TD % 25 m, because in this region there are neither advective nor swimming motions to bring a predator/prey pair into close proximity. The corresponding N; P; Z ð Þprofiles are shown in Fig.  14. The outstanding feature of these results is the formation not of a DBM, but instead a near uniform biological profile extending across nearly all the turbulent depth. But the circumstances of the evolution of this feature are significant. Since N t50 50:1 is so low, the initial behavior of both the P; Z ð Þ fields is to die off through lack of food resources. As the boundary layer is replenished with nutrient (notice it never quite reaches the starting value of N51:0 for the previous simulations), the phytoplankton is able to respond with renewed growth after about 15 d. However, at this time the zooplankton population is still virtually extinct, so there is no predation pressure to retard the phytoplankton growth in the upper regions of the mixing layer. Consequently, the growth is close to uniform across the mixing layer and no significant DBM forms. The decline in the phytoplankton growth rate due to reductions in the ambient light level across the mixing layer is insufficient, in itself, to generate a DBM (There is some variation due to the light. After 18 d, when the global concentration of P reaches a maximum, hPi 51:78 at the surface, increasing to hPi52:06 at z520:65z TD , a difference of 14%. But this is small amount compared to the 50-300% differences seen in the experimental results of Fig. 1 and the simulations in which predation pressure is significant). Instead the strength of the vertical mixing ensures that all the different P cells are exposed to roughly the same amount light, and consequently they all grow at the roughly the same rate. After 20 d or so, the zooplankton population starts to recover in response to the new food resource now available. It too grows uniformly across the boundary layer. One would anticipate that if the simulation were to be run for a longer period of time, then on the next population cycle a situation more akin to the previous simulations would pertain and DBM formation would occur. This is because the nutrient conditions would no longer be so extreme as to force both planktonic fields to near extinction and hence the higher predation pressure near the surface Þforcing. An inward nutrient flux boundary condition is implemented at the surface z50 m of the boundary layer. Biological data as given in Table 2, incorporating the small zooplankton predator, except here it is nonswimming r Z 50:0 m s 21 and fast reacting T R 51 s.
would again manifest itself. It also suggests that the observational DBM recorded by M13 formed as they did because (1) their boundary layers were either well stocked, or recently replenished, with significant amounts of nutrient and (2) the low fluorescence measurements found in the surface regions are a consequence of high predation pressure in those localities. To investigate these points further, Fig. 15 shows what happens in the scenario when the predation term is switched off altogether. Starting from a standard, nutrient rich N; P; Z ð Þ t50 5 1:0; 0:5; 0:5 ð Þinitial profile, the zooplankton population starves and quickly dies out. Unrestrained by predation, the phytoplankton population grows rapidly to very high concentrations, almost uniformly across those parts of the boundary layer where there is sufficient nutrient to In (A), the initial nutrient profile forms a nutricline, ranging from 0:2 to 1:4 for z5233 m to z50 m, respectively. In (B), this initial profile is inverted, with low concentration at the surface and high at the base. Other data as in Fig. 14. support the excess growth (in the top layers for the surface flux, Fig. 15A; at or below z TD for the base flux, Fig. 15B). Again there is little sign that the decline in the ambient light level with depth is a sufficient stimulus (by itself) to initiate significant DBM formation in these wind driven boundary layers.
A variation on this theme is to initiate a simulation with a non-uniform N t50 profile. Two such examples are shown in Fig. 16. Figure 16A shows a run beginning with a moderately high N51:4 nutrient concentration near the surface, declining uniformly to a minimum value of N50:2 at z5z S . The subsequent behavior of the N; P; Z ð Þ profiles from this initial state is interesting. Near the surface the abundance of nutrient promotes rapid growth in the phytoplankton population, but this in turn is rapidly consumed by the zooplankton. Since the reaction time T R 51 s of this small predator is rapid, it is able to feed very efficiently near the surface. Consequently, the phytoplankton population has hardly any time to establish itself before it is consumed and hence the P concentration falls away. One can infer its transient presence by the small predator characteristic s oscil % 6:5 d period still to be seen in the upper zooplankton profiles. Lower down one sees the formation of a very intense DBM at a depth of about 20 m, somewhat deeper than those shown in Figs. 8, 10, 12. Formation at this depth is brought about by the juxtaposition of the relatively high nutrient levels from the initial profile, remaining in situ due to the lack of mixing, in combination with a rapid decline in predation pressure to near zero, again caused by the lack of any mixing or swimming to bring about predator/prey contacts. This DBM also lasts longer than those generated in other simulations and is not subject to the s oscil % 6:5 d period which pertains in the upper layers where strong predation pressure is a factor. The initially imposed nutricline in effect creates two somewhat disconnected boundary layers, in which the biological dynamics evolves at different rates. Figure 16B shows what Þdistribution of (A) the vertical velocity field w x; t ð Þ 5w 00 x; t ð Þ; f since hwi50}, (B) nitrate, (C) phytoplankton, (D) zooplankton and at a depth of z510 m. The data was recorded 4.7 d into the intermediate wind U Ã ; U S ð Þ 5 3:0; 3:3 ð Þ simulation, just at the point when the primary DBM is reaching its maximum extent (see Fig. 10B for the corresponding vertical profiles). Here, the sampling depth lies somewhat above the DBM level.
happens if the initial nutrient profile is inverted from highlow to low-high. As in Fig. 15B, the low level of predation at depths z < z TD combined with easily availability of nutrients (as distinct from Fig. 16A), initially produces very high phytoplankton concentrations. But in this scenario, unlike in Fig.  15B, there is sufficient predation pressure to regulate this growth after just a few days. Higher up (and somewhat overshadowed by the growth near z S ) a DBM again develops around 20 m as in Fig. 16A. The higher ambient nutrient concentration in this case means it is actually twice as intense (cf. the contrasting P concentration scales in Fig. 16) as the DBM which develops from the high-low nutricline.
The intensity scales of both DBM shown in Fig. 16 would indicate that they are somewhat artificial constructs (none of the observed DBM recorded by M13 show such a high level of intensity), brought about by an extremely favorable correlation between the initial nutrient profiles and the rapid decline in both mixing and predation levels below z TD =2.
Nevertheless, they serve to highlight the hypothesis that those regions of the water column which are most favorable to DBM formation are those that exhibit (1), sufficiently high nutrient and light levels to promote growth in combination with (2), a significant fall off in the ambient predation pressure. Since sharp reductions in e z ð Þ bring about a corresponding easing in predation pressure at or just below z ML % z TD =2, it is not surprising that simulated DBM are seen to develop around this depth. It would also explain the strong observational correlations, manifest in M13's datasets, between the location of real DBM and those regions of the water column where the energy dissipation rate is in decline (Fig. 1).

Horizontal distribution
Since the LES-NPZ model carries out 3D simulations, it is worth trying to discern any significant patterns in the evolution of the biological fields when viewed across a horizontal Þsimulation, just at the point when the primary DBM is reaching its maximum extent (see Fig. 10B for the corresponding vertical profiles). Here, the sampling depth corresponds to the center of the DBM.
plane. Periodically during the course of each simulation, instantaneous snapshots of the horizontal distribution of the biological fields are recorded at various different depths. These can then be processed to create "movies" of the biological fields as they evolve during a simulation (copies of these "movies" are available from the authors on request). Some examples of these snapshots, taken from the U Ã ; U S ð Þ5 3:0; 3:3 ð Þ surface nutrient flux simulation (q.v. Fig. 10B for the corresponding vertical profiles), are shown Figs. 17-19. Figure 17 is taken at a depth of 10 m after 4.7 d, while Figs. 18, 19 are for a depth of 20.4 m after 4.7 d and 5.9 d, respectively. The significance of these latter two selections is that they represent samples taken from the primary DBM, which from Fig. 10B is seen to reach its maximum intensity at $ 20 m after around 4 d. But first, it is interesting to consider why the DBM does not form at a shallower depth, at say 10 m for example. Figure 17 shows the instantaneous vertical velocity field w at the same moment as the corresponding biological correlations. The w field exhibits the characteristic Langmuir turbulence pattern of up and downwelling zones (McWilliams et al. 1997;Bees et al. 1998;Lewis 2005;Polten and Belcher 2007;Teixeira and Belcher 2010), in this case, rotated clockwise ($ 40 ) from the wind direction (the x-axis) by a combination of the Coriolis forcing and the Stokes drift term (see Lewis and Belcher 2004). This "streakiness" pattern of elongated up and downwelling zones is replicated (but to a lesser extent) in the biological fields. It is particularly noticeable in the distribution of the P field, which exhibits a correlation between high phytoplankton concentrations and vertical downwelling. Since nutrient is being added at the surface in this instance, the downwelling regions act, in effect, as richer local food environs for the phytoplankton, providing the spur they need for significant extra growth (approximately up to 7% more than the horizontal mean concentration). Notice that this all happens relatively quickly, because the N field shows a significant absence of nutrient in the downwelling zones, presumably because the excess has already been utilised to boost the growth of the local phytoplankton population. Assuming an average downwelling velocity of w % 2310 23 m s 21 around Fig. 19. Key as Fig. 18, except on this occasion the data was recorded after 5.9 d. At this time, the primary DBM is coming to an end of its lifespan, as raised levels of zooplankton predation signal its destruction. this depth, it would take about 500 s for any excess nutrient to be advected a distance of 1 m. This would be time enough to produce differentials in the P concentration field of between approximately 1% and 2% (assuming a light reduced P growth rate of $ 3:3310 25 s 21 ), roughly what is observed. The much slower growing zooplankton (l max Z 52:9310 26 s 21 ) are not sensitive enough to respond to stimuli over as (relatively) short a time period as 500 s, and hence the distribution of the Z field is almost completely uniform in the horizontal (cf. the small color-scale changes in Fig. 17D).
Looking at the corresponding horizontal distributions recorded at the exact same point in time, but at a depth of 20 m (Fig. 18), one is immediately struck by the fact that the vertical mixing is somewhat weaker and the "streakiness" pattern of elongated up and downwelling zones is less pronounced. As a result those regions of the P field distribution exhibiting higher than average growth (shown in yellow), now cover a broader area than in Fig. 17. What seems to be happening is that by the time one reaches a depth of z % 20 m, the excess nutrient is no longer being confined to the strictly defined downwelling zones (Polten et al. 2005) that transported it from the surface. This idea is supported by the velocity variance data shown in Fig. 3A (from a simulation with comparable wind forcing) which shows that at z % 10 m the vertical hw 2 i T velocity component is dominant, but by z % 20 m the horizontal combination of hu 00 2 i T 1 hv 00 2 i T exceeds that of hw 2 i T . This would produce a pooling effect, caused by the replenishment nutrient spreading out (horizontally) from the disintegrating downwelling zones. This pooling of nutrient, in combination with the reduced zooplankton predation pressure at this depth, then promotes the formation of the primary DBM in this region. Notice that the pooling of nutrient is not directly apparent from the N-field itself because on arrival it would be quickly absorbed by the phytoplankton based on the timescales mentioned above. Rather, it manifests itself in the relatively uniform excess growth seen in the P field distribution. In future, it would be interesting to test this pooling hypothesis quantitatively, by carrying out simulations to monitor the average horizontal flux of nutrients at different depths to see if DBM formation corresponds to increased levels of horizontal mixing. Figure 19 shows values of the biological fields at this depth after about 6 d, just as the primary DBM is starting to break down. Over the interval of 1.3 d between Figs. 18, 19, the zooplankton has had the opportunity to generate excess growth in response to the increased higher P field concentrations, increasing the local predation pressure, which ultimately destroys the DBM feature.

Conclusions
The study of the formation and the underlying factors that drive DCM=DBM formation has a long history (Cullen 2015 and reference therein). The literature would seem pretty much exhaustive. However, the recent publication by M13 of various open ocean datasets which unequivocally link DCM formation and levels of background turbulent mixing (something previously that has only been speculated on), provides a new slant and calls for a more detailed investigation of this phenomena. Advances in computing power and resources allow one, using the LES methodology, to simulate 3D, wind and wave driven turbulent boundary layers in great detail. By coupling such simulations to a generic type of NPZ model specially adapted to reflect the influence of background turbulence levels on growth and predation rates, one has the means to carry out just such an investigation.
The results reported here largely corroborate the findings of M13. Provided the wind forcing is not too strong (q.v. Figs. 8C, 10C) and the boundary layer is relatively nutrient rich, the simulations predict DCM/DBM formation at depths at or just below half the turbulent depth % z TD =2, the level at which the upper mixing layer starts to peter out, very similar to the observational data. This depth is a robust feature, as the DCM/DBM continue to form here irrespective of whether the nutrient profile is uniform (primary DBM), or develops in the presence of a nutricline (either bottom up or top down) as seen in the secondary DBM (Figs. 10A,B , 12, 13, 16). The simulations also support the idea that DCM/ DBM are generated primarily in response to predation pressure (Figs. 14-16). Starting the simulations from a very low nutrient base (see Fig. 14) effectively switches off the predation pressure by starving the zooplankton population to near extinction. In its absence, when phytoplankton growth was re-initiated by replenishment of nutrient through the surface, growth resumed uniformly and no DBM formed. Removing the predation term altogether from the model, also inhibited DBM formation (Fig. 15). Taken together, these results indicate that while the vertical mixing is sufficiently vigorous to circulate phytoplankton cells throughout the upper mixing layer so that they all get to experience pretty much the same amount of light, it is not vigorous enough to offset the retardation in growth experienced near the surface through increased levels of predation pressure (see the various predation rate profiles). This is not too surprising because increases in turbulent mixing generate corresponding increases in predation pressure (provided the predator can react fast enough to capture the extra prey contacts it makes), but has relatively little influence on the amount of light penetrating the mixing layer. It is also significant that the DBM formed irrespective of the type of predator (large or small) was used. This conclusion concerning the importance of predation pressure in DCM/DBM formation was also reached by Fennel and Boss (2003) in their (somewhat simpler) mathematical modeling work. It is just possible that DCM/DBM formation could be initiated by a combination of retarded light level growth and a strong positive nutricline (low nutrient at the surface, high at the base) without the need for any predation pressure (this hypothesis is not specifically tested here). However, given the results shown in Fig. 14, such a combination of stimuli does not seem sufficient to produce the intense DBM observed in Fig.  1. Rather, this work tends to support the ideas expressed by Banse (2013) of the importance of "top-down" processes in DCM/DBM formation, since phytoplankton growth is usually matched or exceeded by grazing losses.
The question as to why the DBM form at the depth they do is intriguing. From the investigations here, the depth z TD =2 marks the point at which the coherence of the Langmuir cells driven by the wind and wave generated turbulence starts to break down. The vertical velocity component, which up that point is the dominant feature of the boundary layer turbulence, weakens (see Fig. 3A) and the horizontal components become more prevalent. This leads to a rapid increase in the vertical nutrient gradient as hw 2 i T declines (see Fig. 10), and a smoothing out of the horizontal gradients as the relative strength hu 00 2 i T 1hv 00 2 i T increases (cf. Figs. 16, 17). In combination, this would tend to lead to a pooling of nutrient at around this depth, initiating a strong phytoplankton growth in response. If so, this work indicates that DCM/DBM formation can actually be influenced by the surface wave characteristics driving the Langmuir turbulence regime. A surprising connection, given that in the past it has been assumed that surface wave effects were simply too small to significantly influence even the ocean boundary layer dynamics, let alone the associated biology.
This effect could also be brought about by other features of the flow, such as the presence of shallow thermo/pycnocline acting to restrict the extent of the vertical transport. Since the boundary layers generated in these simulations are unstratified, it is difficult to definitively rule this out. However, in their observations, M13 found that the position of the DCM depth lay some $ 18 m (on average) below the corresponding e z ð Þ peak depth when stratification was relatively low ("A brief . . . Mac ıas et al." section). This strongly suggests that their DCM were not brought about by the presence of a pycnocline. Indeed, Fennel and Boss (2003) argue that in contrast to relatively small lakes (e.g., Simmonds et al. 2015), the density changes associated with shallow pycnoclines found in the open ocean would be too small to disrupt vertical transport very significantly. This idea is supported by the experimental e z ð Þ profiles of M13 (Fig. 1) and the numerical vertical velocity profiles shown in Fig. 9, which both show that vertical mixing does not shut off abruptly at around z TD =2 but continues to exert an influence to a least twice this depth. This would tend to disrupt and deepen pycnoclinic formation, making it unlikely this is the cause of M13's observed DCM.
The DCM/DBM generated in these simulations are transient, lasting no more than a few days. However, this is due to the fact that the biological P2Z fields are two generic representations of many different planktonic species. In reality, DCM/DBM found in the ocean will contain many different species each subject to their own specific growth and predation cycles. So in the absence of any very strong mixing events, observed DCM are likely to be longer lasting than those generated here, although their species composition may change over time. This idea could be tested in more detail by introducing more separate biological species (e.g., by having the large predator feeding directly on the small predator which in turns feeds on the P field) into the simulations, provided suitable biological coupling terms could be introduced into the model equations (Eqs. 3-5). Future developments of this kind should enable better parameterizations for the simpler, but very much faster 1D biophysical models most commonly in use. And if experimental development is such that the recording of physical measurements (e.g., more e z ð Þ sampling) in conjunction with biological ones becomes the norm, then there is scope to develop specialized LES-NPZ models in order to study very specific planktonic ecosystems.