Analytical solution of nitracline with the evolution of subsurface chlorophyll maximum in stratified water columns

In a stratified water column, the nitracline is a layer where the nitrate concentration increases below the nutrient-depleted upper layer, exhibiting a strong vertical gradient 20 in the euphotic zone. The subsurface chlorophyll maximum layer (SCML) forms near the bottom of euphotic zone, acting as a trap to diminish the upward nutrient supply. Depth and steepness of the nitracline are important measurable parameters related to the vertical transport of nitrate into the euphotic zone. The correlation between the SCML and the nitracline has been widely reported in the literature, but the analytic 25 solution for the relationship between them is not well established. By incorporating a piecewise function for the approximate Gaussian vertical profile of chlorophyll, we derive analytical solutions of a specified nutrient-phytoplankton model. The model is well suited to explain basic dependencies between a nitracline and a SCML. The

Abstract. In a stratified water column, the nitracline is a layer where the nitrate concentration increases below the nutrient-depleted upper layer, exhibiting a strong vertical gradient in the euphotic zone. The subsurface chlorophyll maximum layer (SCML) forms near the bottom of the euphotic zone, acting as a trap to diminish the upward nutrient supply. Depth and steepness of the nitracline are important measurable parameters related to the vertical transport of nitrate into the euphotic zone. The correlation between the SCML and the nitracline has been widely reported in the literature, but the analytic solution for the relationship between them is not well established. By incorporating a piecewise function for the approximate Gaussian vertical profile of chlorophyll, we derive analytical solutions of a specified nutrient-phytoplankton model. The model is well suited to explain basic dependencies between a nitracline and an SCML. The analytical solution shows that the nitracline depth is deeper than the depth of the SCML, shoaling with an increase in the light attenuation coefficient and with a decrease in surface light intensity. The inverse proportional relationship between the light level at the nitracline depth and the maximum rate of new primary production is derived. Analytic solutions also show that a thinner SCML corresponds to a steeper nitracline. The nitracline steepness is positively related to the light attenuation coefficient but independent of surface light intensity. The derived equations of the nitracline in relation to the SCML provide further insight into the important role of the nitracline in marine pelagic ecosystems.
The nitracline depth physically depends on the degree of water-column stratification and the magnitude of momentum transfer associated with wind stress (Denman and Gargett, 1983;Laanemets et al., 2004). It also depends on momentum transfer from below (Lipschultz et al., 2002) and, in some cases, vertical advection such as upwelling (Laanemets et al., 2004). However, in a relatively stable environment, the SCML may restrict the diffusive flux of nitrates to the euphotic zone and continually erode the nitracline supposing that sufficient light is available (Probyn et al., 1995). The SCML thereby acts as an effective nutrient trap, regulating the nitracline depth (Banse, 1982;Beckmann and Hense, 2007;Klausmeier and Litchman, 2001;Probyn et al., 1995). However, variation of the nitracline steepness, which is critical to determine the nitrate supply, was poorly understood due to a lack of high vertical resolution data; e.g., both bottle data and Argo data tend to have low vertical resolution sampling. Some studies showed that nitrogen flux is dependent more on the nitracline steepness than on the density gradients regulating turbulent diffusion (Bahamón and Cruzado, 2003;Bahamón et al., 2003;Lavigne et al., 2015). Thus, these measurable features of the nitracline and their correlation with the SCML may provide insightful information for mechanisms of the productivity in pelagic ecosystem, and the analytic solutions for these parameters may fill the knowledge gap.
Although a close relationship between the nitracline and the SCML is always observed, the quantitative nature of the nitracline in relation to the SCML formation has not been studied. The system of phytoplankton and the limiting nutrient on the vertical axis has often been utilized to study the depth, intensity and persistence of the SCML. Major theoretical results include photoacclimation (increase in chlorophyll per cell) (Steele, 1964;Fennel and Boss, 2003), bistability (Yoshiyama and Nakajima, 2002;Ryabov et al., 2010), oscillating SCM (Huisman et al., 2006), hysteresis conditions (Kiefer and Kremer, 1981;Navarro and Ruiz, 2013) and the ESS (evolutionary stable strategy) depth obtained by gametheory approach (Klausmeier and Litchman, 2001;Mellard et al., 2011). Recent mathematical studies solved the persistence and uniqueness of the steady-state solution (Du and Hsu, 2010;Hsu and Yuan, 2010;Du and Mei, 2011) and gave rigorous proofs for the abovementioned ESS depth and the game-theory approach (Du and Hsu, 2008a, b). Additionally, several modeling studies have been conducted to quantitatively assess the importance of different physical-biological processes leading to the SCML (Jamart et al., 1977(Jamart et al., , 1979Varela et al., 1994;Klausmeier and Litchman, 2001;Hodges and Rudnick, 2004;Beckmann and Hense, 2007).
Among the studies using the nutrient-phytoplankton model, Klausmeier and Litchman (2001) first analytically derived the vertical nutrient distribution with the development of the SCML. In this model, the concentration of the limiting nutrient was found to be low and constant above the SCML and linearly increasing with depth below this layer in a poorly mixed water column. Building on that model, Mellard et al. (2011) added stratification and surface nutrient input, which can make phytoplankton grow in both the surface mixed layer and deep layer (SCML) simultaneously. Fennel and Boss (2003) determined that the sum of nutrients and phytoplankton at steady state will increase monotonically below the surface mixed layer until it equals the fixed nutrient concentration. By incorporating a generalized Gaussian function for a vertical chlorophyll profile into the nutrientphytoplankton dynamic equation, Gong et al. (2015) found that the steady-state nitrate concentration increased from the upper community compensation depth to the SCML depth. None of the studies, however, focused on the quantitative nature of the nitracline in relation to the SCML in the stratified waters.
In this paper, we modified the nutrient-phytoplankton model by Gong et al. (2015) to study the roles of SCM in reshaping the nitracline. Two additional terms -atmospheric input, which promotes the growth of phytoplankton in the surface mixed layer, and phytoplankton self-shading, which regulates the light penetration -were introduced into the previous model. Accordingly, a piecewise function comprising a constant value within the surface mixed layer and a Gaussian function below this layer was used as a fit to the steadystate vertical chlorophyll profiles simulated by the nutrientphytoplankton model. By incorporating the piecewise function into the nutrient-phytoplankton model, we derived the analytic solutions for the properties of the nitracline and the SCML in steady state, and the relationship between them was examined in response to light availability, surface nutrient input and vertical diffusivity.
2 Definitions and models 2.1 Models

Dynamic equations
We consider the following equations for phytoplankton and nutrient dynamics in stratified waters (Eqs. 1-2), where light and nitrogen are two limiting factors for phytoplankton growth (Fig. 1). The change in phytoplankton at depth z is the balance of the growth and death and the passive moving (sinking and mixing) (Eq. 1). An eddy diffusion coefficient K v redistributes phytoplankton in the water column. Depth z increases toward the seabed. (Blue solid line is the vertical profile of nitrate gradient; green solid line is Chl a concentration as a function of depth; horizontal green dashed lines indicate the locations of the upper and lower SCML, z m − σ and z m + σ ; horizontal black dotted line indicates the depth of the surface mixed layer, z s ; vertical dotted black line represents zero nitrate gradient; z n1 and z n are the locations of extreme nitrate gradients, z n is also the nitracline depth, and z m is the depth of maximum chlorophyll concentration.) where P denotes the chlorophyll concentration (mg m −3 ). We assume the chlorophyll distribution to represent the distribution of phytoplankton biomass (that means that the photoacclimation of phytoplankton is ignored, and the SCM refers to the subsurface biomass maximum, SBM). This is a significant simplification. In fact, phytoplankton increases intercellular pigment concentration when the light level decreases (Cullen, 1982(Cullen, , 2015Fennel and Boss, 2003). Usually, the depths of the SCML and the SBM layer (SBML) are separate, and the latter is shallower than the former. Nitrogen (in units of mmol N m −3 ) taken up by phytoplankton includes three sources, i.e., recycling from dead phytoplankton, atmospheric input to the surface mixed layer and supply by mixing from deep water (Eq. 2). γ is the nitrogen content of the phytoplankton (mmol N per mg Chl). Following Mellard et al. (2011), we consider the nitrogen input from the atmosphere at the rate of N in (z), setting it as a delta function at z = 0 with the total nutrient input in the surface mixed layer N inML ; N inML = z s 0 N in (z)dz, where z s is the depth of the surface mixed layer. Note that nitrogen input through the activity of nitrogen fixers is excluded. However, Trichodesmium, if it is at the surface, will be modeled similarly to the atmospheric input term.
µ m is the maximum growth rate of phytoplankton, ε is the loss rate of phytoplankton (including respiration, mortality, zooplankton grazing) and α is the recycling efficiency of dead phytoplankton (0 < α < 1). The specific rate of loss processes (ε) is assumed to be vertically homogeneous due to a lack of data (similar to Sverdrup, 1953). The treatment of grazing loss, is, at the least, an oversimplification, though many numerical models used a similar treatment (e.g., Klausmeier and Litchman, 2001;Fennel and Boss, 2003;Huisman et al., 2006). The growth-limited function for light I and nitrogen N is given by µ m min(f (I ), g(N)). The Michaelis-Menten form for the light-limiting function f (I ) and nitrogen-limiting function g(N), f (I ) = I /(K I + I ) and g(N) = N/(K N + N) is used, where K I and K N denote the half-saturation constants for light and nitrogen, respectively. The net growth rate, µ m min(f (I ), g(N))−ε, is positive only if both the light-limiting term µ m f (I ) and nitrogen-limiting term µ m g(N) are larger than the loss rate ε.
Light intensity I decreases exponentially with depth according to the Lambert-Beer law: where I 0 is the surface light intensity and K w and K c are light attenuation coefficients of water and phytoplankton, respectively. For the sake of simplicity, we assume that both K w and K c are constant with depth. The sinking velocity of phytoplankton w is nonnegative in the chosen coordinate system. We assumed it to be constant with depth, neglecting the influencing of density gradients (pycnoclines), which may cause vertical variations in sinking.
To describe the water-column stratification, we assume that the vertical eddy diffusivity K v depends on depth: where z b is set to the bottom boundary of this model and is assumed to be sufficiently deep where the chlorophyll concentration approaches 0. We assume that surface diffusivity (K v1 ) and subsurface diffusivity (K v2 ) (Lande and Wood, 1987;Hodges and Rudnick, 2004) are constant and K v1 is large enough to homogenize the chlorophyll and nitrogen  Raven and Richardson (1986 concentration in the surface mixed layer. A gradual transition from the surface mixed layer to the deep one (in the vicinity of z s ) can be written in terms of a generalized Fermi function (Ryabov et al., 2010), K v (z) = K v2 + K v1 −K v2 1+e (z−zs)/ l , where parameter l characterizes the width of the transition layer. In our numerical study, we assumed this transition layer is 2 m. For simplicity, parameter l is set to infinitely thin in analytic solutions. A comprehensive list of symbols is given in Table 1.
The zero-flux boundary condition for the phytoplankton at the surface is used. At the bottom boundary of the model domain (z = z b ), the Dirichlet or fixed concentration boundary condition is used, i.e., P →0 for z → z b . Fennel and Boss (2003) used an infinite depth as z b (z b → ∞). For the nitrogen distribution we set N = N inML at the surface and diffused into the water column a prescribed nitrate gradient, n, at the bottom. That is,

Fitting equation for vertical chlorophyll profile
In many stratified water columns, the vertical distribution of chlorophyll concentration is homogeneous within the surface mixed layer and appears as a Gaussian below this layer ( Fig. 2a), which is typical in open oceans (Uitz et al., 2006), shelf seas (Sharples et al., 2001), stratified estuary (Lund-Hansen, 2011) and Arctic waters (Martin et al., 2012). The nonuniform vertical profile of chlorophyll within an SCML was first modeled by a generalized Gaussian function (Lewis et al., 1983), which has subsequently been widely used with small modifications. For example, Platt et al. (1988) superimposed a constant background on the generalized Gaussian and fitted it to field data on the vertical distribution of chlorophyll from coastal, upwelling, open oceans and Arctic waters. Afterward, some studies introduced a parameter to represent the slope of the Gaussian curve (Matsumura and Shiomoto, 1993;Mu Oz Anderson et al., 2015). In particular, to account for the observed characteristic that surface values always exceed the bottom ones ( Fig. 2a), the generalized Gaussian functional form has been modified with a superimposition of a background linearly or exponentially decreasing with depth (Uitz et al., 2006;Mignot et al., 2011;Ardyna et al., 2013).
For simplicity, to analytically study the role of the SCML in shaping the nitracline, we therefore propose a piecewise function comprising a constant value in the surface mixed layer and below a general Gaussian function (Eq. 6) to approximate the vertical profile of chlorophyll concentration in where P is chlorophyll concentration as a function of depth z, P 0 is the chlorophyll concentration within the surface mixed layer, and P max = h/ σ √ 2π represents the maximum value of chlorophyll below the surface mixed layer. Considering the influence of the surface mixed layer on the chlorophyll vertical distribution, h is less than the total chlorophyll concentration integrated through the water column. Note that the vertical distribution of chlorophyll is an incomplete general Gaussian function below the surface mixed layer (see Fig. 2a). The three Gaussian parameters (P max , z m , σ ) can characterize the SCM phenomenon. Thus, z m is the depth of the maximum chlorophyll (the peak of the bell shape), and σ is the standard deviation of the Gaussian function, which controls the width of the SCML. The upper and lower boundary of the SCML can be defined as z m − σ and z m + σ , respectively, which are at the depths where there is a balance between phytoplankton growth and losses, and thus they reflect the physical-biological ecosystem dynamics associated with the SCML (Beckmann and Hense, 2007;Gong et al., 2015). To explore the SCML in stratified waters, we assume the surface mixed layer is shallower than the upper boundary of the SCML, i.e., z s < z m −σ . Examples of the piecewise function that fitted vertical chlorophyll profiles in the northern South China Sea (SCS) reasonably well can be found in Gong et al. (2014).
The piecewise function approximation (Eq. 6) was evaluated and justified through numerical simulation of the nutrient-phytoplankton system (Eqs. 1-2), which is solved with a semi-implicit time stepping scheme. The vertical resolution is uniform (2 m), extending down to 200 m. We assumed a small uniformly distributed concentration of phytoplankton (P (z, 0) = 0.1 mg m −3 ) and nitrogen (N(z, 0) = 0.1 mmol N m −3 ) as the initial conditions and run the model until it converges to a steady state (The modeling results show that the steady state has no relationship with the initial values of phytoplankton and nitrogen.) We use the biologically reasonable parameter values given in Table 1 to represent the system at Station SEATS (South East Asia Time Series) in the northern SCS. Thus, the specific (calibrated) model solution is considered as an example to obtain the analytic solutions of the nitracline. Figure 3 shows the numerically simulated equilibrium distributions of nitrogen, light and chlorophyll. In addition, the simulated vertical profile of chlorophyll is fitted well by the piecewise function of chlorophyll using the least square method (Fig. 3). Many numerical solutions of the nutrientphytoplankton system have reproduced the vertical chlorophyll profile with the SCML (Fennel and Boss, 2003;Huisman et al., 2006;Ryabov et al., 2010). Thus, analogous  to the study by Klausmeier and Litchman (2001), we incorporate the piecewise function (Eq. 6) into the nutrientphytoplankton system (Eqs. 1-2) at steady state to examine the roles of the SCML in reshaping the nitracline. We note that the useful delta function approximation in Klausmeier and Litchman (2001) was verified by both simulation and rigorous mathematics (Du and Hsu, 2008a, b). As presented above, the assumption of the piecewise function approximation is physically practical.

Definition of the nitracline
The vertical distributions of nitrate often exhibit a strong gradient in depth (the nitracline), but the features of the nitracline (depth, steepness) are variable in euphotic zones due to the combined effect of physical and biological processes. Many studies define the nitracline depth as the location where the maximum vertical gradient in nitrate concentrations occurs (Eppley et al., 1979;Bahamón et al., 2003;Wong et al., 2007;Beckmann and Hense, 2007;Martin et al., 2010). To measure the defined depth, a high vertical resolution of nitrate concentrations is needed, and this has been a big technical challenge for a long time. Thus, some definitions were also proposed to make the depth measurable. For example, one definition is the depth where the nitrate concentration reaches a prescribed concentration, e.g., 0.05, 0.1, 1.0, or 12 mmol N m −3 (Cullen and Eppley, 1981;Koeve et al., 1993;Martin and Pondaven, 2003). Some studies choose it to be the first depth where nitrogen is detectable (e.g., 0.05 or 0.1 mmol N m −3 ) (Cermeno et al., 2008;Hickman et al., 2012) or where the nitrogen concentration exceeds the mixed layer value by a prescribed concentration difference (e.g., 0.05 mmol N m −3 ) (Laanemets et al., 2004). Significant differences exist between these defined depths, i.e., the depth of the maximal nitrate gradient was found to be deeper by 10 m from the first depth where nitrate can be detected (Eppley et al., 1978), while the nitrate gradient at the first detectable depth of nitrate is nearly 0 (Cermeno et al., 2008).
With the development of nearly continuous nitrate profile measurement using the in situ ultraviolet spectrophotometer (ISUS) optical nitrate sensor (Johnson and Coletti, 2002;Johnson et al., 2010), the detection of the maximum nitrate gradient could be more accurate than before. In this study, we adopt the location of the maximum nitrate gradient (max(dN/dz)) in the euphotic zone as the nitracline depth (z n ), which can be expressed by d 2 N dz 2 z n = 0 and d 3 N dz 3 z n < 0. Below the surface mixed layer, the steady-state version of Eq. (2) reduces to γ µ m min(f (I ), g(N))P − γ αεP = K v2 d 2 N dz 2 . Thus, according to our model approach (Eq. 2), the nitracline depth where d 2 N dz 2 = 0 represents a balance between the nutrient uptake and the recycling of phytoplankton loss, i.e., µ m min(f (I ), g(N)) = αε.
In this study the nitracline steepness is defined as the nitrate gradient at the nitracline depth ( dN dz z n ) (Laanemets et al., 2004;Aksnes et al., 2007).

Data sources
The nitrate profiles were obtained from the ISUS measurement at the SEATS station during the CHOICE-C (Carbon Cycling in China Seas -budget, controls and ocean acidification) 2012 summer cruise. Nine casts were conducted during 6-7 August 2012. The raw ISUS nitrate data, which employed temperature compensation, were first calibrated by the AutoAnalyzer 3 (AA3) and then smoothed to remove noise. The sampling frequency was set to 5 Hz, and the raw data were thus smoothed with a 25-point moving average in the surface mixed layer, a 5-point moving average in the SCML and a 15-point moving average below the SCML. The data were then interpolated by a cubic spline function. The corresponding temperature was obtained from conductivitytemperature-depth (CTD) measurements. Overall, nine sets of profiles are available to examine our analytical solutions.  1) and (2) leads to By substituting the expression of eddy diffusivity (Eq. 4) and the fitted, depth-dependent function of chlorophyll (P (z), Eq. 6) into Eq. (7), we obtain the diffusive term of nitrate below the surface mixed layer, that is, Letting d 2 N/dz 2 = 0 in Eq. (8), for P > 0 one gets Solving this quadratic equation of depth z, we obtain the depths z n1 and z n2 : Taking the derivative of Eq. (8) with respect to depth z, we get K v2 Obviously, at depth z n1 , d 3 N/dz 3 > 0, and at depth z n2 , d 3 N/dz 3 < 0. That is, z n2 is the location of maximum nitrate gradients. We find that the nitracline depth refers to the depth z n2 , i.e., In particular, Eq. (9) becomes a linear function of depth z when the second-order item coefficient (K v2 /σ 4 ) is 0, and thus the equation has only one solution. In fact, in typical stratified waters the diffusivity below the surface mixed layer (K v2 ) is 1-9 × 10 −5 m 2 s −1 , and the thickness of the SCML (2σ ) is from several meters to tens of meters (Cullen, 2015); thus, K v2 /σ 4 (values from 8.64 × 10 −9 to 7.78 m −2 s −1 ) can be neglected in some cases. When K v2 /σ 4 → 0, for non-zero sinking velocity, we get one solution from Eq. (9): Both Eqs. (11) and (12) show that the nitracline depth is located below the SCML depth, i.e., z n > z m (Fig. 1). A numerical study in weak vertical mixing environments showed a similar result (Beckmann and Hense, 2007). Note that the SCML represents the SBML in our model. In some oligotrophic oceans, the SCML will be deeper than the SBML due to the effect of photoacclimation on the vertical distribution of chlorophyll (Fennel and Boss, 2003). For example, Li et al. (2015) showed that the modeled maximum nitrate gradient (nitracline) occurred below the depth of the SCML in the northern SCS, and then we can deduce that the nitracline depth is also deeper than the depth of the SBML. In the Mediterranean Sea, Bahamón et al. (2003) found that the nitracline occurred below the depth of the SCML at 88 % of the stations (50 out of 57 stations). As is well known, the SCML in the Mediterranean is often due to photoacclimation. It is not surprising that the other 12 % of the stations (7 out of 57) in the Mediterranean Sea contradict the finding that z n > z m because photoacclimation leads to a much deeper SCML than the SBML. Thus, we conclude that the nitracline depth is located below the SBML depth.

Nitracline steepness and the SCML
To illustrate the relationship between the nitracline steepness and the SCML, by integrating Eq. (8) from depth z n to z b and using the assumption for phytoplankton at the bottom boundary, i.e., P →0 for z → z b (Eq. 5), we get This equality indicates that the nitracline gets steeper as the distance between the depth of the nitracline layer and that of the SCML is increased. Incorporating Eq. (11) into Eq. (13) leads to Equation (14) indicates that the nitracline steepness is negatively related to the thickness of the SCML.

Depth of the nitracline
By substituting the general Gaussian function for chlorophyll below the surface mixed layer (Eq. 6) into Eq. (1), we obtain the steady-state net growth rate of phytoplankton below the surface mixed layer: Let z 0 = z m − wσ 2 /(2K v2 ) in the first term of the right-hand side of Eq. (15). From the result given by Gong et al. (2015), we know that z 0 is the location of the maximum growth rate of phytoplankton (hereafter referred to as the depth of optimal growth), where the transition from nutrient limitation to light limitation occurs (i.e., f (I ) = g(N) at depth z 0 ). Clearly, z n > z 0 (Eq. 11). Hence, the growth of phytoplankton at the nitracline depth z n is limited by light; i.e., µ m min(f (I ), g(N)) z n = µ m f (I ) z n . In other words, the growth rate of phytoplankton at the nitracline depth is a function of the light level at the nitracline depth, I (z n ). Thus, from Eqs. (11) and (15), we obtain the growth rate of phytoplankton at the nitracline depth, that is, Note that the derivation of Eq. (16) Rearranging Eq. (17), we find µ m − αε = αεK I /I (z n ). This equality indicates that the maximum rate of NPP, (µ m − αε)P , is inversely proportional to the light level at the nitracline depth, I (z n ). Lande et al. (1989) found that higher maximum rates of population growth corresponded to shallower nitracline depths in the central North Atlantic. Furthermore, the insertion of Eq. (3) into Eq. (17) yields another expression of the nitracline depth: Note that Eq. (18) is obtained on the premise that the nitracline depth exists. This equality shows that the nitracline depth is inversely proportional to the light attenuation coefficient of water (K w ), and it deepens logarithmically with increasing surface light intensity (I 0 ). It is noted that the nitracline depth has a negative relation with the self-shading of phytoplankton (K c z n 0 P dz). Importantly, Eq. (18) predicts that the nitracline depth has no relation with subsurface diffusivity. Aksnes et al. (2007) also proposed a similar result, i.e., that a shoaling nitracline per se cannot be taken as an unequivocal sign of increased upwelling, as well as eddy diffusion. However, this does not mean that fluid dynamics are unimportant in shaping the vertical distribution of nitrate.
Equation (18) also indicates that both a higher recycling efficiency (α) and a larger loss rate (ε) lead to a shallower nitracline, while the enhanced maximum growth rate of phytoplankton (µ m ) moves the nitracline depth downward. Modeling results showed that the nitracline was shoaled by 24 % (from 84 m upwards to 64 m) when both the recycling efficiency (α, from 0.6 to 0.8) and the loss rate (ε, from 0.3 to 0.4 d −1 ) were increased by 33 %. Accordingly, the predicted nitracline depth from Eq. (18) varied from 86 to 71 m. Increasing µ m by 33 % (from 0.9 to 1.2 d −1 ) makes the simulated nitracline deepen slightly (from 84 to 88 m), leading to the predicted nitracline depth changing from 86 to 92 m. The experiments with varying parameter values indicate that the updated z n (based on the model runs) matches the predicted z n of Eq. (18) well.

Steepness of the nitracline
In steady state, integrating Eq. (2) from the nitracline depth z n to the bottom boundary z b , and considering the light limitation of phytoplankton growth below depth z n (Eq. 15), we obtain the nitrate gradient below the surface mixed layer: This equality shows that the nitracline steepness enhances with increasing nitrate gradient at the bottom boundary ( dN dz | z b ) which depends on the intensity of nitrate intrusion from below. The vertical diffusion negatively influences the nitracline steepness. The modeled time series distributions of nitrate gradients and diffusive nitrate fluxes in the northern SCS and the upstream Kuroshio Current showed similar results (Li et al., 2015). Beckmann and Hense (2007) conducted a sensitivity analysis of both vertical diffusivity and nutrient concentration at the bottom boundary to examine the vertical phytoplankton and nutrient profiles in oligotrophic oceans, and their numerical results support the relations presented in Eq. (19).

Analytical solutions for SCM characteristics
Similar to methods used by Gong et al. (2015), the piecewise function for the vertical chlorophyll profile (Eq. 6) was incorporated into the nutrient-phytoplankton model (Eqs. 1-2) at steady state to derive the three SCM characteristics (SCML thickness, its depth and intensity).
For z = z m and z = z m + σ , the net growth rate of phytoplankton (Eq. 15) can, respectively, be expressed as follows: By substituting the growth limitation function for light to Eq. (20) or Eq. (21), we obtain the expression of parameter z m , i.e., Subtracting Eqs. (22) and (23), and rearranging, we obtain the expression of parameter σ : Neglecting terms including self-shading of phytoplankton (K c ) in Eqs. (22-24), both the analytical solutions of the depth and thickness of the SCML are the same as the results presented in Gong et al. (2015). The self-shading effect of phytoplankton plays an important role in the vertical pattern of chlorophyll (Shigesada and Okubo, 1981;Klausmeier and Litchman, 2001;Beckmann and Hense, 2007). In line with common sense, our analytic results indicate that a higher selfshading of phytoplankton negatively influences the depth and thickness of the SCML (Eqs. 22 and 24), having a similar effect as an increasing light attenuation coefficient of water, K w . The expression of the SCML intensity is different from the results presented in Gong et al. (2015). Integrating Eq. (7) from the surface of the water to the bottom of the surface mixed layer (z s ) and from the bottom of the surface mixed layer to the base of our model domain (z b ) gives: Parameter λh (mg m −2 ) is assumed as the integrated chlorophyll concentration below the surface mixed layer. Based on the properties of the Gaussian function, λ can be expressed as λ = can be obtained from the standard normal table. According to the property of the Gaussian function, we have 0.68 < λ < 1.00, z s < z < z b , under the assumption of z s < z m − σ .
Adding Eqs. (25) and (26) yields Nitrogen input to the surface mixed layer (N inML ) causes an increase in surface chlorophyll concentration (Eq. 25). Hence, the total chlorophyll in stratified water columns (λh+ P 0 z s ) increases with increasing N inML (Eq. 27), which has also been predicted by the numerical study of Mellard et al. (2011) and supported by the experimental test of Mellard et al. (2012). Because recycling processes are assumed to not immediately convert dead phytoplankton back into dissolved nutrients below the surface mixed layer, i.e., α = 1, the total chlorophyll concentration below the surface mixed layer and the intensity of the SCML can, respectively, be expressed as: The integrated chlorophyll concentration below the surface mixed layer (λh) and the intensity of the SCML (P max ) are influenced positively by N inML and negatively by P 0 (Eqs. 28-29). That is to say, the influence of nitrate input to the surface mixed layer on the SCML intensity (also on the integrated chlorophyll concentration below the surface mixed layer) is nonlinear. Hence, their changes (λh and P max ) with varying N inML cannot be obtained from the steady-state solutions straightforwardly, depending on the specific parameter combinations in the model. For example, λh and P max decrease when increasing nutrient enrichment directly to the surface mixed layer in the ecosystem given by Mellard et al. (2011), while they are nearly unchanged in oligotrophic oceans (Varela et al., 1994). Our results (Eqs. 28-29) also show that enhanced subsurface diffusivity (K v2 ) increases the integrated chlorophyll concentration and the intensity of the SCML (λh and P max ), as a result of a higher nitrate flux (K v2 n). Physical upward transport of nitrate across the bottom of the nitracline is indeed the main nitrogen source for NPP in the euphotic zone (Ward et al., 1989).

Light effects on the nitracline with the SCML
We now examine how the steady-state nitracline, in relation to the SCML, depends on light availability, especially the light level at the nitracline depth.
Substituting Eq. (17) to Eq. (28) and rearranging, we have This equality indicates that the light level at the nitracline depth, I (z n ), is positively related to the integrated chloro-phyll concentration in the whole water column, λh+P 0 z s . we can find from Eq. (3) that the nitracline depth (z n ) is inversely related to integrated chlorophyll. This inverse relationship has been observed in many regions. In southern California coastal waters, the phytoplankton standing stock and its primary production rate were positively related to the reciprocal nitracline depth (Eppley et al., 1978(Eppley et al., , 1979. Bahamón et al. (2003) found that larger depth-integrated chlorophyll with an average deeper SCML and nitracline (∼ 129 and ∼ 136 m, respectively) occurred in the western Sargasso Sea, central Sargasso Sea and eastern Atlantic than in the Canary Current zone.
The nitracline depth deepens with increasing surface light intensity but with decreasing light attenuation coefficients (K w and K c ). These results were consistent with observations; e.g., Letelier et al. (2004) found the depth of the nitracline to coincide with an isolume, a depth of constant light level in the North Pacific Subtropical Gyre.
The predicted effect of surface light intensity and light attenuation coefficient on the nitracline depth (Eq. 18) implies that the nitracline depth in stratified waters may have seasonal variations. In the North Pacific Subtropical Gyre, Letelier et al. (2004) found that the nitracline depth differences between winter and summer disappeared when nitrate concentrations were plotted against light level in the water column. Aksnes et al. (2007) found that the seasonal pattern of nitracline depth was governed by seasonality in the light attenuation coefficient, rather than in surface light intensity. In particular, the inverse proportional relationship between the nitracline depth and light attenuation coefficient (Eq. 18) has also been derived from a steady-state model by Aksnes et al. (2007), which is consistent with observations in the coastal upwelling region off Southern California (Aksnes et al., 2007). Teira et al. (2005) found a significant positive correlation between the nitracline depth and the depth of 1 % surface light intensity (the proportion of reciprocal light attenuation coefficient) in the eastern North Atlantic Subtropical Gyre. Bahamón et al. (2003) showed that the nitracline depth remained relatively constant around 1 % surface light intensity depth in the western Sargasso Sea.
The nitracline steepens with a higher light attenuation coefficient (K w and K c ) due to K w and K c negatively influencing SCML thickness (Eqs. 14 and 25). Numerical modeling showed that a higher K w leads to a thinner SCML and thus a steeper nitracline layer (Beckmann and Hense, 2007). Aksnes et al. (2007) also found that the fluctuations in the nitracline steepness were positively correlated with the fluctuations in reciprocal Secchi depth (i.e., light attenuation coefficient) in the upwelling area off the coast of southern California. We further point out that the nitracline steepness almost stays constant with changing surface light intensity (I 0 ) because surface light intensity has no relation to the SCML thickness (Eq. 25). The sensitivity analysis of a onedimensional (vertical) model showed that the vertical nutri-ent profiles were almost parallel to each other with increasing surface light intensity (Beckmann and Hense, 2007).
The inverse effects of the light attenuation coefficient on the nitracline steepness and its depth imply that the nitracline becomes steeper as it shoals. Aksnes et al. (2007) found this consistent pattern in the upwelling area off the coast of southern California.

In the presence of surface nutrient input
Current evidence and modeling analyses suggest that climate warming will increase ocean stratification and hence reduce nutrient exchange between the ocean interior and the upper mixed layer (Cermeno et al., 2008;Chavez et al., 2011). Therefore, nutrient input directly to the euphotic layer due to atmospheric deposition may become a relatively more important nutrient supply mechanism to the euphotic layer (Mackey et al., 2010;Okin et al., 2011;Mellard et al., 2011). However, few model studies (e.g., Mellard et al., 2011) have explored the influence of external surface nutrient supply on vertical phytoplankton distribution.
Observations show that an interzone exists between the transition of the surface mixed layer and the deep layer, where the nutrient gradient equals nearly 0 dN dz z s +0 = 0 (Fig. 3), leading to the solution in Eq. (26). It follows that the total chlorophyll in the surface mixed layer depends on the surface nutrient supply (N inML ) (Eq. 26). In this case, if N inML is negligible, Eq. (26) degenerates to (1 − α) εγ P 0 z s = 0.
In this case, the dead phytoplankton in surface mixed layer must be fully recycled, i.e., α = 1, in order to sustain the positive chlorophyll concentration (P 0 > 0). In other words, if the dead phytoplankton cannot be fully recycled in the surface mixed layer, an external nutrient supply to the layer is needed to fuel the growth of phytoplankton. Thus, the term external nutrient supply to the surface mixed layer should be included in the system equations at steady state to make a positive surface chlorophyll concentration. Numerical results by Mellard et al. (2011) also showed that phytoplankton populations can grow in the mixed layer and in the deep layer together, when there is nutrient input directly to the mixed layer. However, a surface nutrient source is not a necessary term for a model approach where dissolved organic matter and detritus are explicitly resolved (Beckmann and Hense, 2007). Accordingly, we treat the vertical phytoplankton distribution as a piecewise function, comprised a linear function in the surface mixed layer and a Gaussian function below, which is more realistic than the general Gaussian function.
The assumption of the piecewise function for phytoplankton is also consistent with the assumption of piecewise vertical diffusivity. For simplicity, we assume that the transition layer between the surface mixed layer and the deep one is infinitely thin and that the chlorophyll is continuous within the transition layer. By assuming that the SCML depth is significantly deeper than the base of the surface mixed layer, we obtain the steady-state solutions for the SCML depth and thickness, similar to the solutions using the general Gaussian function. However, the intensity of the SCML is affected by a surface nutrient supply with an associated positive increase in phytoplankton concentration.

SCML trapping nutrients
Observations and numerical simulations showed that the SCML played a role as a nutrient trap in some regions, restricting the diffusive flux of nitrates to the surface mixed layer (Anderson, 1969;Klausmeier and Litchman, 2001;Navarro and Ruiz, 2013).
From Eq. (10), we know z s < z n1 < z 0 −σ < z m −σ (Fig. 1). That is, the SCML occurred below depth z n1 . For z m < z n (Eq. 10), we know that the upward diffusive nitrate is sufficient to sustain the fast growth of phytoplankton in the lower part of the SCML (z m < z < z m +σ ). To explore the SCML restricting nitrates into the surface mixed layer, we next examine if the nitrate concentration at depth z n1 above the upper boundary of the SCML (z n1 < z m − σ ) is depleted.
According to the definition of the depth z 0 (where f (I ) = g(N) holds) and z 0 > z n1 (Eq. 10, Fig. 1), we know that the growth of phytoplankton at depth z n1 is nitrate-limited, i.e., µ m min(f (I ), g(N)) z n1 = µ m g(N) z n1 . From Eq. (14), we find that at depth z n1 , the growth rate of phytoplankton equals the recycling efficiency of dead phytoplankton, i.e., µ m g(N) z n1 = αε. Inserting the Michaelis-Menten form for g(N) into this equality yields N(z n1 ) = K N /(µ m /αε − 1). Phytoplankton maximum growth rates (µ m ) of 0.2 to 1 per day are typical in optical environmental conditions (Banse, 1982;Timmermans et al., 2005). We choose 0.5 per day to illustrate the result. Loss rate (ε), although not well documented, is often quoted as 10 % of the growth rate (Parsons et al., 1984). A reasonable choice for the remineralization efficiency seems to be α = 0.5 (Huisman et al., 2006). The typical value of half-saturation constants for nitrate (K N ) is between 0.1 and 0.7 mmol N m −3 in oceans (Eppley et al., 1969). We adopt 0.4 mmol N m −3 . Thus, we find that the nitrate concentration at depth z n1 , N(z n1 ), is equal to 0.03 mmol N m −3 , a value lower than the detection limit, indicating the depletion of nitrate above depth z n1 .
Because the SCML acts as a nutrient barrier, it is easy to understand that the rate of NPP in the SCML ((µ m min(f (I ), g(N)) − αε)P , z m − σ < z < z m + σ ) is positively related to upward nitrate flux that is trapped. This condition can simply be derived by integrating Eq. (2) vertically at steady state, i.e., z m +σ z m −σ (µ m min(f (I ), g(N)) − αε)γ P dz = K v2 dN/dz z m +σ z m −σ . This result suggests that, the production within the SCML is fueled mainly by nitrate and is thus NPP. Because, at the nitracline depth, the gross growth rate µ m min(f (I ), g(N)) equals the recycling of dead phytoplankton αε, we assumed for the constant αε that within the nitracline layer (z n1 < z m −σ and z m < z n2 ), the nitrate uptake by phytoplankton has to be supplied by the vertical diffusion. Observations also showed that most of the primary production in the SCML was supported by nitrate from vertical diffusion, with an average f ratio (i.e., relative contribution of the nitrate uptake to the total nitrogen uptake) of 0.74 ± 0.26 during early summer in Canadian Arctic waters (Martin et al., 2012).

Vertical profile of nitrate gradients
From the monotonicity of the quadratic function of depth z in the left-hand side of Eq. (9), we know that d 2 N/dz 2 < 0 when z s < z < z n1 and z n2 < z < z b , but d 2 N/dz 2 > 0 when z n1 < z < z n2 . In other words, the vertical gradient of the nitrate concentration (dN/dz) decreases with depth on the interval (z s , z n1 ), while it increases on the interval (z n1 , z n2 ) and then decreases on the interval (z n2 , z b ). If we consider the distribution of vertical nitrate gradients as continuous across the base of the surface mixed layer, then we get dN/dz < 0 for z s < z < z n1 under the assumptions of the uniform nitrate distribution within the surface mixed layer (i.e., dN/dz = 0, 0 < z < z s ). The schematic of vertical profiles of nitrate gradients and chlorophyll concentrations in stratified waters is shown in Fig. 1.
The negative gradient of nitrate below the surface mixed layer (dN/dz < 0 for z s < z < z n1 ) indicates that the nitrate concentration decreases with depth on the interval (z s , z n1 ). This decreasing nitrate feature has rarely been observed by traditional measurements probably due to the techniquelimited low resolution. Some float data showed this feature in vertical nitrate profiles; for example, Sakamoto et al. (2009) found it at depths below the base of the surface mixed layer (∼ 45-50 m) using the ISUS temperature-compensated data at an eastern Pacific oligotrophic station. Our in situ time series measurements using the ISUS at SEATS station also showed this decreasing feature at depths ∼ 25-30 m (Fig. 4). We note that this decreasing nitrate feature will disappear in our derivation if the subsurface vertical diffusion is too weak (Eq. 12) or the surface mixed layer is deeper than depth z n1 . Simulating results showed that the negative gradient of nitrate becomes smaller with increasing sinking velocity (w) and recycling efficiency (α). This implies that the negative gradient is likely the result of the consumption of nitrate by phytoplankton exceeding the supply of nitrate.

Limitation and application
The model in this study integrates a number of physical, chemical and biological processes that act together to determine the vertical distribution of phytoplankton and nitrate, under the assumption that the system is strictly vertical and in steady state. A few processes such as oxygen status, photoacclimation, luxury uptake of nutrients, phytoplankton motility, concentration-dependent herbivory and depth-dependent herbivory are not included, although they can affect the vertical distribution of phytoplankton and nitrate. Detritus, dissolved organic matter and zooplankton are not included explicitly, and all loss processes, except for sinking, are set to be linearly proportional to phytoplankton. The sinking velocity of phytoplankton is assumed to be independent of density gradients. Further, the vertical transport of nutrients is only by eddy diffusion in our model; in reality, nutrients can be supplied by many processes (turbulence, internal waves, storms, slant-wise and vertical convection), especially by upwelling (Katsumi and Hitomi, 2003;Aksnes et al., 2007).
In this study, the sinking velocity of phytoplankton is set to be independent of nitrate concentration. A vertically varying sinking velocity has been observed as a physiological response to variations in light or nutrient levels (Steele and Yentsch, 1960;Bienfang and Harrison, 1984;Richardson and Cullen, 1995). The sinking velocity reduces with decreased light level and with increased nutrient concentration, and the resulting divergence in sinking velocity can be large enough to affect the location of the phytoplankton particle maximum. However, numerical results given by Fennel and Boss (2003) showed analytically that the divergence of the sinking rate contributes to the location of the SBM layer in a significant way only when the divergence in sinking rates occurred above the compensation depth in stable, oligotrophic environments. They also found that in stable, oligotrophic environments with a predominance of small cells, the biomass maximum is located at the depth where growth and losses are equal, leaving little influence by sinking divergence.
It is worth pointing out that, in extreme oligotrophic regions, the SCML is very deep and attributable mostly to photoacclimation of chlorophyll content rather than to a peak in biomass (Steele, 1964;Fennel and Boss, 2003;Cullen, 2015). The process of photoacclimation is also important for the nutrient-phytoplankton system (with stratified conditions) we focused on. To explore the influence of photoacclimation on the nitracline, we parameterized Chl : C using the mathematical description by Cloern et al. (1995), i.e., Chl : C = 0.003+0.0154 exp(0.050T ) exp(−0.059I )µ . That is, Chl : C is the ratio of Chl a to C in phytoplankton growth at steady state under defined temperature T ( • C), daily irradiance I (mol quanta m −2 d −1 ) and nutrient-limited growth rate µ (µ = N/(K N + N)). The ratio Chl : C increases when temperature T or nitrogen concentration N increases, while it decreases with increasing daily irradiance I . Let R = Chl : C; then, the nitrogen content of phytoplankton can be written as γ = 1/(6.625 × 12 × R), corresponding to a C : N ratio of 6.625 and a carbon atomic mass of 12. From the expression of nitracline depth (Eq. 18), we know that the ratio Chl : C has no influence on the nitracline depth, while the nitracline steepness increases with increasing parameter γ (Eq. 19). In other words, the nitracline gets steeper with a lower ratio of Chl : C. Note that a model with equations that explicitly resolve variations of the Chl a-to-carbon and nitrogen-tocarbon ratio of phytoplankton would certainly be more realistic.
The piecewise equation (Eq. 6) can be used to mimic a large variety of vertical chlorophyll profiles from coastal, upwelling, open oceans and high-latitude waters (Fig. 2). For example, for z s > 0, when the depth of the surface mixed layer equals or is deeper than the depth of the SCML, the vertical profiles like Fig. 2b and c are often found in well-mixed waters (Uitz et al., 2006). For z s = 0, the vertical distribution of chlorophyll concentration (Fig. 2d) can be expressed by a Gaussian function, which is usually found in coastal upwelling waters (Xiu et al., 2008). In particular, when z s = 0 and z m = 0, the surface bloom occurs (Fig. 2e). In general, the vertical profiles of chlorophyll can be classified as belonging to two types, i.e., one peak distribution and one uniform distribution in large regions of lakes and oceans (Uitz et al., 2006;Lavigne et al., 2015). Note that the skewed profiles of chlorophyll with a sharp SCM were not considered in this study. The small-scale (1 m or less) vertical heterogeneity in chlorophyll distribution has been shown to be common features in coastal waters (Sullivan et al., 2010;Prairie et al., 2011;Durham and Stocker, 2012), referred to as a thin layer.
Choosing the values of model parameters representative of the system in the northern South China Sea (given in Table 1), we can retrieve the nitracline depth and steepness, the optimal depth, and the three SCM characteristics. To make calculation easy, we neglect the term of self-shading by phytoplankton in the calculation because a higher self-shading parameter has the same effect as an increasing light attenuation coefficient by water. The calculated and observed values of these parameters are listed in Table 2. All these parameters calculated are in a reasonable range, although there are some discrepancies compared with observations. In fact, this is not surprising, considering that we assume a single phytoplankton group and neglect the microbial loop and the dynamics of the dissolved organic matter and detritus pools. We stress that the analytical solutions of the nitracline are valid only for estimates of z m , h and σ that are consistent with the model's numerical steady-state solution. The numerical steady-state solution in turn depends on the combination of parameter values and on the forcing and boundary conditions. The approximations of z m , h and σ are entirely conditioned by the modeling results and thus also depend on the combination of model parameter values. To combine the analytical steady-state solutions with observed z m , h and σ (as derived from vertical profiles of chlorophyll a concentration) is only meaningful after model calibration (identifying a model solution that is in some agreement with the observed z m , h and σ ).

Summary
We have presented a theoretical framework to investigate the interaction of phytoplankton and nutrient in a stratified water column. A piecewise function for chlorophyll profiles comprising a linear function in the surface mixed layer and a Gaussian function below is assumed in the nutrientphytoplankton model at steady state. A number of important findings are obtained under the conditions of the model equations imposed.
In steady state, the nitracline is confined to between two depths where the gross growth rate equals the recycling efficiency of dead phytoplankton, indicating that within the nitracline, nitrate consumption by phytoplankton has to be replenished by the upward flux of nitrate. This layer thereby is the major contributor to new primary production.
The nitracline depth is located below the SCML depth; both depths deepen with either increased surface light intensity or a decreased light attenuation coefficient. The nitra-cline depth does not depend on the value of the subsurface diffusivity. The nitracline is steeper with a thinner SCML. The nitracline steepness is positively influenced by the light attenuation coefficient, yet it responds insignificantly to surface light intensity.
Our analytical solutions show that phytoplankton in the SCML acts as an efficient nutrient trap, filtering out the upward nitrate supply. The light level at the nitracline depth has a positive relation with the depth-integrated chlorophyll concentration in the whole water column and with the maximum rate of NPP, acting as the indicator of integrated NPP. The NPP is constructed from the model equations that rely on Blackman's law of limiting factor for the growth rate.
Data availability. The data used in this study are available upon request to the corresponding author.
Competing interests. The authors declare that they have no conflict of interest.