The shape of the oceanic nitracline

In most regions of the ocean, nitrate is depleted near the surface by phytoplankton consumption and increases with depth, exhibiting a strong vertical gradient in the pycnocline (here referred to as the nitracline). The vertical supply of nutrients to the surface euphotic zone is influenced by the vertical gradient (slope) of the nitracline and by the vertical separation (depth) of the nitracline from the sunlit surface layer. Hence it is important to understand the shape (slope and curvature) and depth of the oceanic nitracline. By using density coordinates to analyze nitrate profiles from autonomous Autonomous Profiling EXplorer floats with InSitu Ultraviolet Spectrophotometers (APEX-ISUS) and shipbased platforms (World Ocean Atlas – WOA09; Hawaii Ocean Time-series – HOT; Bermuda Atlantic Time-series Study – BATS; and California Cooperative Oceanic Fisheries Investigations – CalCOFI), we are able to eliminate much of the spatial and temporal variability in the profiles and derive robust relationships between nitrate and density. This allows us to characterize the depth, slope and curvature of the nitracline in different regions of the world’s oceans. The analysis reveals distinguishing patterns in the nitracline between subtropical gyres, upwelling regions and subpolar gyres. We propose a one-dimensional, mechanistic model that relates the shape of the nitracline to the relative depths of the surface mixed layer and euphotic layer. Though heuristic, the model accounts for some of the seasonal patterns and regional differences in the nitrate–density relationships seen in the data.


Introduction
Dissolved inorganic nitrate is an essential macronutrient for oceanic phytoplankton production.In the vast majority of the surface ocean that is sunlit, the production of phytoplank-ton is limited by the availability of nitrate.In such regions, consumption of nitrate by phytoplankton renders the sunlit ocean devoid of nitrate.Nitrate increases rapidly with depth in the underlying region, termed the nitracline, where it is resupplied through the remineralization of organic matter and redistributed along isopycnals by the ocean's circulation.The vertical supply of nitrate to the surface ocean depends not only on the vertical transport induced by dynamical processes like turbulent entrainment, Ekman pumping, and frontal and eddy-induced upwelling, but also on the vertical gradient of nitrate.Further, it depends on the depth from which nitrate needs to be drawn.Often, the isopycnal beneath which there is a significant vertical gradient of nitrate, termed the nitrate-depletion density, and its depth, termed the nitracline depth, can be clearly distinguished.This isopycnal is indicative of the density (and depth) from which nitrate must be transported into the euphotic layer, either by (reversible) uplift or by (irreversible) advective and turbulent mechanisms.
One of the fundamental difficulties in finding generalized descriptions for the vertical distribution of nitrate is the high degree of variability in the nitrate-depth profiles.This may be partially overcome by exploiting the nitrate-density relationship that has long been known to be more robust than the relationship between nitrate and depth across the main pycnocline (Redfield, 1944;Pytkowicz and Kester, 1966).In some areas, the nutrient-depleted near-surface layer transitions abruptly to a uniform gradient of nitrate with respect to density (termed the nitracline slope), whereas in others, the transition is more gradual, leading to curvature in the nitratedensity relationship.We term the slope and curvature of this relationship the nitracline shape.The shape and depth of the nitracline are pertinent to the vertical supply of nutrients for new phytoplankton production.The objective of this study is to characterize these properties of the nitracline from data Published by Copernicus Publications on behalf of the European Geosciences Union.and search for patterns and explanations for the revealed relationships.Strickland (1970) was probably the first to find that over weeks, density could be used to predict nitrate to almost within experimental accuracy irrespective of the depth of the nutrient measurement.Recent studies (McGillicuddy et al., 1999;While and Haines, 2010;Ascani et al., 2013) have confirmed that nitrate is better correlated with density than with depth in the ocean.By relating isopycnal excursions with sea surface height, Ascani et al. (2013) used the data from profiling floats to show that a higher variability of nitrate along isobars, as compared to isopycnals, on timescales ranging from hours to weeks, can be ascribed to the movement of isopycnal surfaces by internal waves and eddies.The relationship between nitrate and density also extends to longer time-and space scales, leading to a large-scale alignment between the mean isopycnal and iso-nitrate surfaces (Omand and Mahadevan, 2015).
Since observations of nitrate are relatively rare and more difficult to obtain than measurements of temperature and salinity, a well-characterized relationship between nitrate and density could be useful in a number of applications ranging from the modeling of primary production, where the initialization and restoration of nitrate is crucial, to the interpretation of phytoplankton productivity from space.Several studies have exploited the nutrient-density relationship to develop descriptive algorithms for nitrate.Elrod and Kester (1985) applied a cubic spline method to produce reference curves for nutrient and oxygen in the Sargasso Sea.Kamykowski and Zentara (1986) developed third-order polynomial fits between nutrients and density (and nutrients and temperature) within 10 • regions from the National Oceanographic Data Center (NODC) data set.Garside and Garside (1995) used a stepwise polynomial method to predict nitrate from temperature and salinity during regional phenomena, such as the North Atlantic spring bloom and El Niño/La Niña conditions in the central Pacific.Finally, firstorder linear regressions between nitrate and temperature have been used to predict nitrate in upwelling zones (Traganza et al., 1987;Dugdale et al., 1997;Olivieri and Chavez, 2000;Omand et al., 2012), although mixing and interleaving of water masses may confound this approach (Friederich and Codispoti, 1981).In coastal California, the relationship has been used to predict coastal nitrate concentration over the last century as an indicator of giant-kelp health (Parnell et al., 2010).Algorithms that exploit the nitrate-temperature relationship have also been developed for remotely sensed sea surface temperature (SST) and ocean color (Dugdale et al., 1989;Goes et al., 2000;Switzer et al., 2003); however, the strong correlation between nitrate and temperature tends to break down near the surface as biological effects and upperocean heat fluxes alter the relationship (Garside and Garside, 1995).Biological processes challenge the predictive capacity of the empirical algorithms, and there are relatively few simple mechanistic models for oceanic-nitrate distribution in density coordinates (Kamykowski, 1987;Dugdale et al., 1989;Omand and Mahadevan, 2015).
In this study, we draw on the growing data set of nitrate measurements to obtain a characterization and synthesis of the vertical distribution of nitrate.Vertical profiles of dissolved inorganic nitrate from shipboard sampling during programs such as the Hawaii Ocean Time-series (HOT), the Bermuda Atlantic Time-series Study (BATS), the California Cooperative Oceanic Fisheries Investigations (CalCOFI) and the compilation of climatological gridded data sets such as the World Ocean Atlas (WOA09 Garcia et al., 2010) provide long time series (> 15 years).The more recently developed optical sensor, the in situ ultraviolet spectrophotometer (ISUS) (Johnson and Coletti, 2002), has facilitated the high-resolution in situ measurement of nitrate.By outfitting APEX autonomous profiling floats with such sensors, Johnson et al. (2013) are able to offer detailed sampling of the vertical distribution of nitrate in time (Johnson et al., 2010).Such a growing database of nitrate measurements offers the opportunity to explore the vertical distribution of nitrate in space and time and the regional variability in the shape and depth of the nitracline.Further, we can explore the relationship of the nitrate distribution to other parameters.
Here, we combine empirical and modeling approaches to describe the nitrate-density relationship.We apply secondorder polynomial fits to the above data sets, diagnosing the nitrate depletion density, the slope and an index of curvature for the nitrate-density relationship.We develop an idealized model that predicts the curvature of the nitrate-density profile as a simple one-dimensional function of the euphotic depth and mixed-layer depth.We find that in the subtropical gyres, the shape may be predicted by the interplay between the depth of the mixed layer, the euphotic zone and the remineralization depth.In some regions, however, upwelling, lateral mixing and transport by the large-scale circulation can complicate the relationship, rendering a one-dimensional approach less suitable.Finally, we discuss the applications of this approach to ocean numerical modeling and to predictions of the vertical turbulent nitrate flux at the base of the euphotic zone.

Data sets
We used the publicly available nitrate (NO 3 ) and potential density (σ t ) data from the World Ocean Atlas (WOA09), long-term time series (HOT, BATS, CalCOFI) and APEX-ISUS floats in the Pacific and Atlantic subtropical gyres.Only a brief description of each of these data sets is provided here.
(WOA09) (Garcia et al., 2010), available through the National Oceanographic Data Center (http://www.nodc.noaa.gov/OC5/SELECT/dbsearch/dbsearch.html).The data set compiles roughly 50 years of ship-based profiles of nitrate and density into a climatology by averaging data on a 5 • grid and into 33 depth levels, with a vertical bin size that varies from 10 m near the surface to 400 m below 2000 m depth.

HOT, BATS and CalCOFI
The data set was constructed from the analysis of bottle samples of NO 3 from monthly cruises between 1996 and 2012 at the deep-water station ALOHA (A Long-Term Oligotrophic Habitat Assessment; 22 • 45 N, 158 • 00 W), located 100 km north of Oahu, Hawaii, at BATS Hydrostation(31 • 40 N, 64 • 10 W) in the western Atlantic ocean, and from a 22year record of quarterly measurements compiled from twelve CalCOFI stations (spanning 121.5 to 124.4 • W and 30.5 to 32.7 • N).The CalCOFI profiles are horizontally averaged over the region to generate a single time series of profiles.

APEX-ISUS floats near Hawaii and Bermuda
Nitrate and potential density are obtained from two Webb Research APEX profiling floats (Johnson et al., 2013); one deployed from 4 December 2009 (to present) at the ALOHA station near Hawaii (6401hawaii) and a second, deployed from 6 November 2009 to 1 November 2011 near the BATS station (6391bermuda).Each float collected a profile of conductivity, temperature, pressure, NO 3 and oxygen between the surface and the float's parking depth of 1000 m, with measurements vertically spaced every ∼ 5 m (between 7 and 100 m depth) and 10 m (between 100 and 400 m), at roughly 5-day intervals (www.mbari.org/chemsensor/FloatList.html).Here, we examine only NO 3 and σ t .Nitrate concentration is estimated with an in situ ultraviolet spectrophotometer (ISUS) with a short-term precision of about ±0.1 µM; however, over multiyear deployments, sensor drift introduces an additional uncertainty, resulting in an accuracy of about ±0.4 µM (Johnson et al., 2010(Johnson et al., , 2013)).

Methods
The typical profile of NO 3 vs.density in the ocean is described as follows: within the least-dense, sunlit surface layer, nitrate is depleted due to uptake by phytoplankton.Across the main pycnocline, NO 3 increases with density, often reaching a subsurface maximum.In depth coordinates, this subsurface maximum in nitrate generally occurs between 500 and 1000 m, or intersects the seafloor.We focus on the nitrate-density relationship spanning the region of the main pycnocline that lies above the deep NO 3 maximum and beneath the NO 3 depletion density (Fig. 1a).

Second-order polynomial fitting
The inspection of NO 3 -σ t profiles from various regions over the range described above reveals that while some profiles are very linear, others are not (Fig. 1b-g).We therefore seek a second-order polynomial of the form NO σ 3 = a (σ t −σ o ) 2 + b (σ t − σ o ) + c to describe the data from each region or location and then compare the coefficients to distinguish regional characteristics.We subtract the density (σ o ) at the upper limit of the fit to minimize any effect in the fitting due to the mean value.Although using a higher-order polynomial would improve the fit skill, we restrict ourselves to a secondorder polynomial because our goal is to characterize the most basic features of the NO 3 -σ t relationship in an effort to gain a mechanistic understanding of the processes that govern the nitracline shape and depth.In addition, for a number of the ocean time series and WOA09 profiles, poor vertical resolution means that the inclusion of three or more parameters in the functional fit would increase the risk of overfitting the data.The coefficients a, b and c in the polynomial fit are obtained through a least-squares minimization of where k denotes each level in a profile or areal average of profiles, n is the number of data points in the profile and NO σ 3 denotes the nitrate fit (in contrast to NO 3 , which we use to denote data).The skill of each fit is evaluated over the bounded region, according to the normalized mean rms error calculated as The slope of NO σ 3 (σ t ) is given by dNO σ 3 dσ t = 2a σ t − b.The "linearity" of the profile can be estimated as the difference in the local slopes of NO σ 3 at either end of the profile (2aσ t −b)| k=n k=1 = 2a(σ t,k −σ t,1 ) = 2a σ t , as compared with the slope over the entire profile, The ratio of these gives a dimensionless curvature index (the non-dimensionalized coefficient a) given by where NO σ 3 and σ t denote the change in nitrate and density at either end of the profile.A large value of a arises from a nonlinear profile where the NO 3 gradient varies locally with σ t .A small value of a implies a linear profile with a constant gradient.Hence, we define three categories of fits.Those that are nonlinear (a > 1), those that are linear (a < 1) and those that are poorly fit (skill< 0.8, due to higher-order terms or scatter).
In the circumstance in which the fits are linear (a < 1), we also fit NO σ 3 (σ t ) with a straight line of the form NO   (4) We can then find the nitrate depletion density, σ depl , where NO σ 3 goes to 0 as σ t o = −c lin /b lin + σ o .Similar to Kamykowski and Zentara (1986), we suggest that σ depl represents the deepest isopycnal at which nitrate is depleted (or the nitrate depletion density).Since σ t tends to be more highly resolved than discretely sampled NO 3 , this provides an appealing way of estimating the nitracline depth z o in circumstances where NO 3 is poorly vertically resolved or unavailable.

Fitting considerations
Particularly for nonlinear fits, the coefficients are sensitive to the upper and lower bounds of NO 3 that we select.For each profile, we aim to fit only the region beneath NO 3 depletion and above the NO 3 maximum with positive curvature.However, due to the large number of profiles, we tried to establish general criteria for the upper and lower extents of the fitted profile segments.As a shallow limit on the NO 3 fits, we retain only data with NO 3 > 2 µM.This threshold generally distinguishes the NO 3 -depleted surface samples in the discretely sampled time series and WOA09 data and falls sufficiently above the 0.4 µM accuracy of the ISUS sensors on the floats (Johnson et al., 2013).The depth of this threshold varies spatially (Fig. 2a) and temporally, particularly in loca-tions with strong seasonality in production, upwelling, (e.g., the California Current) or mixed-layer depth (e.g., Bermuda).The deep bound on the NO 3 -σ t fits is qualitatively defined as some σ t above the deep-NO 3 maximum.For the time series data (long-term stations and APEX-ISUS floats), we found that the NO 3 maxima tends to remain associated with a particular water density over time (for example, see Figs. 6 and  7).We select thresholds of σ t < 26.75, 26.5 and 27.2 kg m −3 that fall roughly 0.5 kg m −3 above the NO 3 maximum at Cal-COFI, HOT and BATS, respectively.At HOT, for example, this threshold was selected near the inflection point where NO 3 begins to curve downward (e.g., see Fig. 6).The same thresholds as at HOT and BATS are used for the corresponding APEX-ISUS float profiles.In the WOA09 climatology, the selection of a the deep threshold was more challenging, since each 5 • × 5 • grid has different profile characteristics.Thus, we manually selected the lower bound on the fit for each grid by selecting the inflection point for each NO 3 -σ t profile individually.The inflection point is demonstrated in Fig. 1a, and the results for the lower bound on the WOA09 fits are shown in Fig. 2b.
The number of data points that are retained between the shallow and deep bounds described above can also affect the fit coefficients, particularly a .Through bootstrapping both the climatological and time series data, we found that in general, the sensitivity of the results to a were significantly diminished after the inclusion of five or more points in each fit.In the WOA09 climatology, the upper and lower bounds were selected to include at least five levels (Fig. 2c).The average number of data points used in the time series fits were 14 ± 5 (WOA09), 15 ± 8 (CalCOFI), 7 ± 5 (HOT), 8 ± 4 (BATS), 10 ± 2 (6391bermuda) and 17 ± 3 (6401hawaii).For each fit, we evaluate the skill according to Eq. ( 2).Where the skill is poor (operationally defined here as < 0.8), the profiles are excluded from the analysis.In general, we find that the NO 3σ t skill exceeds 0.8 for nearly all of the WOA09 data set (Fig. 2d) and 90 % of the time series profiles.

Results
In the following sections, we present results from the fitting method to describe global NO 3 and σ t climatology (WOA09 5 • ×5 • ) and time series data from the California Current (Cal-COFI), the North Atlantic Subtropical Gyre (BATS, APEX-ISUS float) and the North Pacific Subtropical Gyre (HOT, APEX-ISUS float).We analyze the fit coefficients (a, b, c), a , b lin and c lin to explore the spatial and temporal patterns and variability in these different records.We develop a simple 1D model for predicting a from a balance between the mixed-layer depth (z ML ) and euphotic depth (z eu ) and compare this with the subtropical APEX-ISUS data.

Nitrate-density fits from the WOA09 climatology
The nutrient depletion density σ depl (Fig. 3a ogy (Fig. 3b-d), is similar in pattern to the nitrate-depletion temperature described in Kamykowski and Zentara (1986) according to cubic regressions between NO 3 and temperature in 10 • gridded climatology.Since this 1986 publication, roughly 700 000 new observations of NO 3 have been added to the National Oceanographic Data Center database.For a comparison with the results of Kamykowski and Zentara (1986), we perform the same analysis between climatological WOA09 NO 3 and temperature (T ) measurements.We find that the skill in NO 3 -T (not shown) has a very similar magnitude and distribution to the σ t fits (Fig. 2d), and our method is not able to distinguish whether T or σ t is a better proxy for NO 3 .We also compute a for each 5 • ×5 • grid (Fig. 4a).In 59 % of these cases, we find that a < 1 indicates a linear NO 3 -σ t relationship between the shallow and deep bounds of our fits.We observe a large-scale spatial coherence of a in some regions.For example, we find that a in the eastern subtropical Pacific, throughout most of the Atlantic and in the Northern Indian Ocean is predominantly less than 1.In the Southern Ocean and the equatorial Pacific (incidentally, regions that are less nitrate limited), a tends to be larger than 1, indicating curvature in NO 3 -σ t .The map of a shows large variation between adjacent grid cells that are not attributable to known biogeographic boundaries.We explore the possibility that these results for a were instead related to the number of points used in the fits (Fig. 2c) or lower skill (Fig. 2d) but these factors do not explain the variability.The WOA09 climatology is compiled from profiles collected at different times of year and from the number of total samples.It is likely that some of the variation between grids is a result of inhomogeneous sampling effort.
Where a < 1, we find the slope (b lin ) according to Eq. (4) (Fig. 4b).In these locations (representing 59% of the climatology), the shape of the NO 3 -σ t relationship between nitrate depletion and the nitrate maximum is σ t -independent and described by a single parameter b lin .The slope varies from roughly 2 µM m 3 kg −1 , up to 18 µM m 3 kg −1 at very high latitude.Some of the low values in the western tropical Pacific, for example, may be attributed to the large change in density (stratification) between the surface and depth of the nitrate maximum (≈ 150 m; Fig. 2b).Conversely, some higher values of b lin correspond to regions where the upper bound on the NO 3 fits are deep (> 150 m; Fig. 2a).

Nitrate-density fits from the California Current region (CalCOFI)
Several previous studies have identified a linear nitratetemperature (NO 3 -T ) relationship beneath the nitratedepletion density in the California Current region (Strickland, 1970;Parnell et al., 2010;Omand et al., 2012).Profiles of NO 3 -σ t in the spatially averaged quarterly CalCOFI record appear relatively consistent over time (green segments, Fig. 5a).We compute a for each profile and find that a is consistently less than 1, suggesting that the NO 3 -σ t relationship is linear over the 22 years analyzed (black points, Fig. 5b), and consistent with the WOA09 analysis, which suggests that the relationship is relatively linear throughout the eastern North Pacific Ocean below 50 • N (Fig. 4a).Interannual variations are dominant over seasonal and spatial variability in a (Fig. 5b).The slope of the linear NO 3 -σ t fit (4) is shown for the annually and spatially averaged record (Fig. 5c).The nutrient depletion density (σ depl ) varies interannually from 24.6 to 25.4 kg m −3 (black circles, Fig. 5d).
The linear slope b lin is correlated with σ depl (r 2 = 0.61, p < 0.001) and is inversely correlated (r 2 = 0.21, p < 0.001) with the nutrient depletion temperature (T depl ; red circles, Fig. 5d).We hypothesize that variations in b lin and σ depl may be related to the SST Nino index (yellow bars, Fig. 5d), which assesses warming in the surface ocean of the eastern Pacific.Indeed, some years do appear to have an aboveaverage T depl (the 1998 El Niño) or low σ depl (the 2004-2007); however, the pattern does not hold for all cases, indicating that there are additional drivers of the NO 3 -σ t relationship.

Nitrate-density fits from HOT and float 6401hawaii
Monthly σ t versus bottle-sampled NO 3 from the HOT station at 22 • 45 N, 158 • 00 W (Fig. 6a) is compared with a higher depth-and time-resolved σ t versus NO 3 record from a nearby APEX-ISUS float (6401hawaii), drifting over a nominally 400 km 2 region, centered at 23 • N, 161 • W (Fig. 6b).The water column in this oligotrophic region remains stratified year-round, with near-surface (nominal 7 m) density fluctuations between 24 and 24.6 kg m −3 , seen both in the HOT and APEX-ISUS records.Nutrient depletion (NO 3 < 2 µM) occurs near the 25 kg m −3 isopycnal, and the nitrate maximum occurs near 27.1 kg m −3 (in 2012) over the 2.8-year record.Profiles of NO 3 -σ t show a marked difference in shape from those in the California Current.Over the fitted sections of each profile, NO σ 3 (σ t ) has positive curvature, with a > 1 for 201 of the 206 total APEX-ISUS profiles and increasing from 1 to 2 (Fig. 6c).
Nitrate (color, Fig. 6d) within the shallow mixed layer (black line, Fig. 6d) is depleted year-round, and the base of the euphotic layer is consistently much deeper (around 180 m as indicated by the gray lines, Fig. 6d).The depth of the mixed layer, z ML , and euphotic layer, z eu , are obtained as described in the Appendix.

Nitrate-density fits from BATS and float 6391bermuda
Monthly σ t versus bottle-sampled NO 3 from the BATS station at 31 • 40 N, 64 • 10 W (Fig. 7a) and a more highly depth-and time-resolved σ t versus NO 3 record from a nearby APEX-ISUS float (6391bermuda), drifting over a nominally 600 km 2 region, centered at 32 • N, 66 • W, are shown in Fig. 7b.The surface (nominal 7 m depth) density of the Sargasso Sea fluctuates seasonally between 26.5 and 23.5 kg m −3 through deep winter convection and stratification in the summer.This 2-year record is punctuated by the passage of mesoscale eddies and wintertime storm events (Lomas et al., 2013, and references therein).In contrast to the strong seasonality and transient events in the density structure, the overall relationship between NO 3 and σ t over our fitted region varies relatively little (colored regions, Fig. 7a,  b).The average linear slope (b lin ) is 32.3 ± 4kg m −3 µM −1 .The deep nitrate maximum is centered at 27.3 kg m −3 , and nitrate depletion occurs near the wintertime mixed layer potential density σ t = 26.5 kg m −3 .This depletion density cor- responds to the density of the 18 • C subtropical mode water.The production and advection of this water mass may alter downstream nutrient delivery (Palter et al., 2005).
The NO 3 -σ t relationship varies seasonally in both BATS and APEX-ISUS records.It is linear between January and July (a < 1, green segments) and curved between August and December (a > 1, yellow segments).Scatter between successive a tends to be small relative to the seasonal variation (Fig. 7c).Fluctuations in the nitracline depth due to deep mixing (Steinberg et al., 2001) and mesoscale eddies (McGillicuddy et al., 1998), though apparent in depth coordinates (Fig. 7d), are not distinguishable in density coordinates (Fig. 7ab).The mixed-layer depth (z ML ; black line, Fig. 7d) is deepest in late winter (max.380 m) and shoals to less than 10 m in the late summer.The euphotic depth (z eu ; gray lines, Fig. 7d) remains between 100 and 200 m, resulting in |z eu | > |z ML | during the summer and |z eu | < |z ML | for a portion of the winter.In the following section, we propose that in nitrate-limited systems, the nitracline shape may be controlled at first order by the depth of the euphotic layer relative to the mixed layer z eu − z ML .

A one-dimensional model of the nitracline shape
Following Lewis et al. (1986), we model the vertical profile of nitrate assuming a balance between the vertical supply of nitrate to the euphotic layer, the uptake of nitrate by phytoplankton in the euphotic layer and the resupply of nitrate by remineralization.This can be expressed as where κ(z) represents the vertical eddy diffusivity and κ(z) ∂NO 3 ∂z the vertical diffusive flux of nitrate, which is set to 0 at the upper and lower boundaries.The uptake of nitrate is modeled as proportional to the initial slope of the photosynthesis-irradiance curve (α, (µM m 2 W −1 )), the fraction of total nitrogen uptake that is supplied vertically in the form of nitrate (the "f ratio", γ NO 3 ), and the light profile, with the incident photosynthetically available irradiance (E o (W m 2 )) attenuated exponentially according to a diffuse attenuation coefficient k (m −1 ; yellow shaded region, Fig. 8).
The euphotic depth is defined where the light attenuates to 1 % of the incident level (E o ) such that z eu = − log(0.01)k.
The NO 3 uptake rate, µ, is comparable to phytoplankton growth rates (∼ 0.1 to 2 d −1 ).Here, we use a low-growth (uptake) rate µ = α γ E o = 0.17 d −1 characteristic of subtropical gyres (Maranon, 2005).The remineralization of organic nitrogen to NO 3 is modeled as the divergence of the vertical flux and taken to be exponential in depth (green shaded region, Fig. 8), with a vertical scale (k −1 rem ) of 250 m, which is typical of observed (Lutz et al., 2002) and modeled (Kwon et al., 2009) flux profiles.The magnitude of the flux R at z = 0 is selected so that total NO 3 is conserved and the verti-cally integrated loss of NO 3 due to uptake is balanced by the vertically integrated NO 3 that is restored between the bottom of the domain, z = −H , and the surface, z = 0.
We model the density as subject to a weak surface heat flux implemented as κ(z) ∂σ t (z) ∂z| z=0 = C −1 p α 0 q, where C p is the specific heat capacity of water (J-kg −1 • C −1 ), α 0 is the thermal expansion coefficient of water ( • C −1 ) and q is the surface heat flux (q = 2 W m −2 ).
The vertical eddy diffusivity κ(z) is prescribed to be depth-dependent and to enhance the vertical mixing over the mixed layer as where κ min and κ max are the lower (deep) and upper (mixedlayer) diffusivities and z w is a vertical length scale for the transition (here z w = 200 m, though the results are not sensitive to this choice).We select κ max = 10 −4 m 2 s −1 and κ min = 10 −5 m 2 s −1 , values that are similar to Lewis et al. (1986).The modeled κ(z) profile is shown in Fig. 8 (pink shaded region).We solve Eqs. ( 5) and ( 7) over a 1000 m (H ) water column.We set σ t to a constant value at z = −H , and prescribe a weak surface heat flux q, but no surface or bottom flux of NO 3 .In steady state, the divergence of the vertical NO 3 flux is balanced by uptake and remineralization.The two parameters that we vary are the mixed-layer depth, z ML , and the euphotic depth, z eu .We find that the curvature of the NO 3 -σ t relationship depends strongly on these two free parameters.
We implement the model with z ML and z eu chosen to be typical of late winter at HOT (CASE 1: z ML = −30, z eu = −210 m) and BATS (CASE 2: z ML = −210, z eu = −80 m) stations.For CASE 1, σ t = 24.9kg m −3 at the surface, with an upper-and lower-layer buoyancy frequency (N 2 ) of 1.3 × 10 −5 and 5 × 10 −5 s −2 , respectively (black line, Fig. 8a).In CASE 2, we maintain the same κ in the upper and lower layers, but the mixed layer is deeper, and σ t = 27.5 kg m −3 at the surface, reflecting the cooler surface water at BATS (black line, Fig. 8b).NO 3 it depleted near the surface in both cases due to uptake, with a subsurface maximum due to restoring (Fig. 8c, d).However, in CASE 2, the uptake of NO 3 occurs only within the mixed layer, and so the large κ that operates on the density also acts to homogenize NO 3 over the same region, and both NO 3 and σ t have a sharp transition near z ML .In contrast, CASE 1 has a deep euphotic layer, and so NO 3 uptake occurs below z ML , resulting in a smooth NO 3 transition approaching the surface.This difference is reflected in the different shapes of the NO 3 -σ t relationship (Fig. 8e, f) and is quantified using polynomial fits for estimating a from Eq. (3) as was done for the data in the previous section.We fit the solution to a quadratic over the range NO 3 > 0.5 µM and σ t < 26.5 and 27.2 kg m −3 following the lower bounds for HOT and BATS, respectively.NO 3 -σ t for CASE 1, with |z ML | < |z eu |, is curved (a = 1.62;Fig. 8e), whereas for CASE 2, with |z ML | > |z eu |, the profile is linear (a = 0.05; Fig. 8f).These curvatures of the profiles bear resemblance to randomly selected APEX-ISUS profiles from the late winter of 2010 (Fig. 8g, h) in the same regions.We also formulated the model with a light-saturating growth curve but found that the steady-state solution remained unchanged.
We explore the dependance of the modeled a on the parameter z ML − z eu .Keeping z ML = 100 m, κ max , κ min , κ w , k rem and the boundary conditions fixed as for CASE 1, we obtain solutions to Eqs. ( 5) and ( 7) over a range of z eu , from −10 to −300 m.We find a from the polynomial fit between NO 3 and σ t and solutions for each new z eu (central black curve, Fig. 9).We find that the curvature index a is a monotonically increasing function of z ML −z eu , with a 1 where z eu is shallower than z ML .Near z ML − z eu = 70 m, a surpasses 1 and increases steeply until leveling off near a = 2.5.
Our previous examples (CASE 1 and 2; Fig. 8) fall near this continuum on either side of the a = 1 threshold.We test the sensitivity of these results to the µ and κ max parameters by varying µ from 0.05 to 0.8 d −1 and κ max from 2 × 10 −5 to 4 × 10 −4 m 2 s −1 .Increasing µ by a factor of 4 shifts the a curve to the right, suggesting that the NO 3 -σ t relationship remains linear for a deeper z eu than the previous cases (gray lines, Fig. 9).Conversely, decreasing µ by a factor of 4 shifts the curve to the left, indicating that a smaller difference between z ML and z eu is required before a > 1.We find the a is relatively insensitive to the choice of κ max in the range 10 −3 to 10 −5 m 2 s −1 (thin black lines, Fig. 9).

Model-data comparison
Our one-dimensional model predicts that when |z eu | < |z ML |, the nitracline shape is linear in density space (a 1) and when z eu is significantly deeper than z ML , the nitracline is curved (a > 1).The model assumes that the divergence of the horizontal NO 3 flux is small compared with the vertical flux divergence, which may not be applicable in settings where lateral water mass intrusions strongly modify the profiles.It also assumes that at depth H = 1000 m, σ t remains fixed and that NO 3 uptake balances remineralization over this depth range.Acknowledging the limitations of these assumptions, we plot the APEX-ISUS a (black circles, Fig. 10) against the corresponding z ML − z eu (see Appendix).
We find that for 6391bermuda, a varies with z ML − z eu (Fig. 10b) but the relationship is scattered and also contains a hysteresis where the transition from a < 1 to a > 1 is delayed compared to the transition from deep to shallow mixed layers (see Fig. 7c, d).For 6401hawaii (Fig. 10b), the model fails to explain the variability in a -reflected as an overall trend towards more positive a (Fig. 6c).The scatter and independence of a relative to our z ML − z eu index suggests that other processes or parameters are likely also important in determining the nitracline shape.However, the model does predict a > 1 at Hawaii and the much larger variability in a at Bermuda, and the results generally fall within the modelpredicted bounds (Fig. 10).Thus, we conclude that a simple model similar to the one presented here may be a useful framework for understanding some of the basic features and differences in the nitracline shape in different regions, particularly where nitrate is limiting, the horizontal gradients are relatively small and variations in the deep nitrate supply are relatively steady.

Nitrate in ocean models
The primary productivity in ecosystem models, and global carbon cycle models, relies strongly on the distribution of nitrate in the model.Data sets of nitrate are invaluable in initializing these models.If nitrate and density are initialized from concurrent measurements or a consistent database, the nitrate-density relationship is implicit.In many instances, however, it would be advantageous to initialize the model's density field and use a nitrate-density relationship for prescribing nitrate.This way, the NO 3 -σ t relationship is maintained even if the model's density is altered.Treating nitrate in density coordinates is also useful if a restoring scheme is used for nutrients at depth.

Implications for estimating turbulent vertical fluxes of NO 3
The turbulent vertical flux of NO 3 (F NO 3 ) can be written as the product of the vertical eddy diffusivity (κ) and the vertical gradient in nitrate (∂NO 3 / ∂z) according to The vertical eddy diffusivity, κ, is difficult to estimate.It may be inferred from tracers or microstructure measurements of the turbulent dissipation rate, (e.g., Dewey and Crawford, 1987).Based on observations of the relationship between shear production and buoyancy flux, κ is approximated by the relation (Osborn, 1980)  where γ is a mixing efficiency (approximately 0.2), the buoyancy frequency is ∂z and is the turbulent kinetic energy dissipation rate.Applying the chain rule to Eq. ( 8), whereby ∂z , and substituting Eq. ( 9) into Eq.( 8), the N 2 dependence cancels, yielding We expect that F NO 3 will be proportional to and the gradient of NO 3 across isopycnals.The turbulent flux of NO 3 near the base of the euphotic zone is of particular interest as the NO 3 flux into the euphotic zone supports new primary production in many oceanic regions.The striking temporal consistency in the NO 3 -σ t relationship from this analysis, suggests that this relation could prove useful for estimating vertical NO 3 fluxes from , which, unlike κ, can be directly measured.

Depth of the nitracline
One way of defining the depth of the nitracline is to use a nitrate cutoff and test for the depth where the nitrate crosses this value.This may lead to noisy results because nitrate can be highly variable within the euphotic layer, and the vertical resolution of the data may be insufficient to identify the cutoff depth.Aksnes et al. (2007) used a novel approach of relating nitrate in the near-surface ocean to the attenuation profile of light by solving a one-dimensional equation in which the uptake of nitrate by phytoplankton is linearly related to light and nitrate itself.Since light decays exponentially with a characteristic attenuation coefficient, the depth and slope of the nitracline can be related to the attenuation characteristics of light, which are easily measured.Here, we suggest that by fitting the nitrate to density, one can identify the nutrient depletion density σ t,o where NO * 3 = 0 and thereby identify the depth of the isopycnal σ t,o as the depth of the nitracline.This gives robust results when the NO 3 -σ t relationship is linear, in which case σ t,o = b lin /c lin .

Conclusions
Having examined nitrate profiles from several sources, the WOA09 gridded data, BATS, HOTS, CalCOFI time series and float records from Bermuda and Hawaii, we find that the NO 3 -σ t relationship can be characterized using a secondorder polynomial to fit the data.The nondimensional curvature a serves as a measure of nonlinearity of the relationship.When linear, a < 1 and the slope of the profile can be identified, as can the nitrate depletion density and nitracline depth.A simple one-dimensional model is able to explain the nitracline shape (a < 1 or a > 1) in terms of the depth of the euphotic layer, z eu , and mixed-layer depth, z ML .When |z eu | − |z ML | > 100 m, a > 1. Conversely, when |z eu | − |z ML | < 100 m, a < 1.Though the model makes the assumption of one-dimensionality, we find that |z eu | − |z ML | is broadly indicative of whether the NO 3 -σ t relationship is linear or curved in the subtropical gyre data.The model relies on the assumption of weak variation in lateral nitrate supply and thus may not be applicable in more complex regions.It is advantageous to examine NO 3 in an isopycnal coordinate system as compared to depth coordinates because variability due to vertical mixing and advection affect NO 3 and σ t in a similar way while generally maintaining the characteristics of the NO 3 -σ t relationship.With the growing database of nitrate from ship-based and autonomous platforms, a more detailed picture will undoubtedly emerge and provide the opportunity to refine our understanding of the underlying processes that govern the nitracline shape and nitrate supply to the surface ocean.

Appendix A
The mixed-layer depth (z ML ) is computed from monthly CTD (conductivity-temperature-depth) profiles as the depth where the density difference is 0.03 kg m −3 from the nearsurface density; it is then smoothed with a 90-day low-pass filter.At HOT, the pycnocline is present year-round and z ML varies between −20 and −80 m on sub-seasonal timescales (Fig. 6d).At BATS, stratification is eroded during the winter as the surface loses heat, and wind and wave-driven turbulence mix the upper water column to depths ranging between −200 and −400 m.During the summer, an intense pycnocline develops and z ML shoals to less than −20 m (Fig. 7d).
The euphotic depth (z eu ) is defined as the level where photosynthetically available radiation (PAR) is 1 % of the surface PAR.Because in situ PAR is not measured at HOT, BATS or on the APEX-ISUS floats, z eu is estimated two ways: (1) from the depth integral of Chl a following Morel (1988) and Morel and Berthon (1989), henceforward referred to as MB (solid gray line, Figs.6d and 7d), and (2) from the Moderate Resolution Imaging Spectroradiometer (MODIS) level-3 K490 satellite product (dashed gray line, Figs.6d and 7d).The appropriate depth over which to integrate Ch a in the above expression is determined iteratively, first using Ch a integrated between the surface and 400 m.The z eu MB is used as the next limit of integration for Ch a .This process is repeated until the solution for z eu MB converges.
Satellite-derived euphotic depth (z eu sat ) was determined from z eu sat = −4.6/K490, where the 8-day, level-3, MODIS K490 is spatially averaged over a 10 • × 10 • region centered on the BATS and HOT stations.At BATS, both z eu MB and z eu sat vary seasonally, becoming shallower during the winter and deepening during the summer.The z eu MB is based on in situ measurements of Ch a throughout the water column, whereas K490 is a function of the water-leaving radiance over the first optical depth (K490 ≈ z eu /4.6).Both are approximations of the true z eu , and differences between z eu sat and z eu MB reflect the methodological differences.The z eu sat integrates water-leaving irradiance over a shallow portion of the water column, roughly 5 to 20 m thick, whereas z eu MB integrates over the entire water column.
For example, in the spring of 2010 at BATS, a subsurface maximum in Ch a developed that led to shoaling of z eu MB (solid gray line, Fig. 7d).These events were not evident in z eu sat , possibly because this subsurface layer was not visible within the satellite penetration depth.For the model-data comparison in Sect.5.1 (Fig. 11), we use z eu sat .The variation in z ML -z eu at Bermuda is driven primarily by variation in z ML , and so our analysis was not sensitive to this choice.
Figure 1.(a)Schematic of a typical oceanic NO 3 versus σ t profile.In nitrate-limited regions, near-surface (low-density) NO 3 is near 0 above a nitrate depletion density and tends to increase monotonically with increasing σ t until reaching a subsurface NO 3 maximum.The region below NO 3 depletion and above the inflection point (where NO 3 begins to curve downward towards the maximum, black segment) is evaluated with a second-order polynomial (NO σ 3 , yellow line), where a is an index of curvature.Panels (b-g) show examples from the WOA09 climatology (circles), with black circles indicating those that are used in the fits.Colored lines indicate a < 1 (green), a > 1 (yellow) or skill < 0.8 (red).

b
lin (σ t − σ o ) + c lin , where b lin and c lin are obtained from the least-squares minimization of NO 3k = b lin σ k + c lin , k = 1, . ..n.

Figure 2 .
Figure 2. Global maps from the 5 • WOA09 climatology indicating the (a) depth of the upper limit (NO 3 > 2 µM), (b) subsurface NO 3 maximum, (c) number of depth levels that fell between the depths in (a) and (b) and were used in each fit and (d) skill of the polynomial fit to the data evaluated as skill = 1 − 1 n

Figure 4 .
Figure 4. Global maps from the 5 • WOA09 climatology of (a) a and (b) b lin , where a < 1. Green regions (a < 1) in (a) indicate a linear NO 3 -σ t fit.
Figure 5. (a) Successive profiles of NO 3 as a function of σ t from quarterly CalCOFI cruises.Green segments indicate the portion of each profile that is fit with Eq. (1), and profiles without green segments have a fit skill < 0.8.(b) a (black dots) for each of the profiles from (a).Annual averages of the (c) slope b lin and (d) intercepts σ t o (black) and T o (red) of linear σ t -NO 3 fits Eq. (4), with error bars indicating the standard deviation from the spatial and temporal averaging.Yellow shaded regions in (d) indicate the time periods with a positive SST anomaly in the Nino-3.4region associated with El Niño conditions.

Figure 6 .
Figure 6.Successive profiles of NO 3 as a function of σ t from (a) the monthly HOT time series and (b) an APEX-ISUS float (6401hawaii) deployed nearby from December 2009 to June 2012, profiling at 5-day intervals.Green and yellow segments indicate a < 1 and a > 1, respectively (skill > 0.8).(c) a for each of the APEX-ISUS profiles in (b).(d) Depth-resolved time series of NO 3 (colors) with σ t contours overlaid.The mixed-layer depth, z ML , and the euphotic depth, z eu , are shown with the black and gray lines, respectively.The solid and dashed gray lines represent different approximations of z eu (see Appendix).

Figure 7 .
Figure 7. Successive profiles of NO 3 as a function of σ t from (a) the monthly BATS time series and (b) an APEX-ISUS float (6391bermuda) deployed nearby from November 2009 to November 2011 profiling at 5-day intervals.Green and yellow segments indicate a < 1 and a > 1, respectively (skill > 0.8).(c) a for each of the APEX-ISUS profiles in (b).(d) Depth-resolved time series of NO 3 (colors) with σ t contours overlaid.The mixed-layer depth, z ML , and the euphotic depth, z eu , are shown with the black and gray lines, respectively.The solid and dashed gray lines represent different approximations of z eu (see Appendix).

Figure 8 .
Figure 8. Examples of profiles of σ t (a and b) and NO 3 (c and d) produced with the simple one-dimensional model.The mixed-layer depth (z ML ) and the euphotic depth (z eu ) are chosen to represent typical conditions near Hawaii (with a shallow z ML = −30 and deep z eu = −210 m; upper panels, "CASE 1") and near Bermuda during late winter (with a deep z ML = 180 and shallow z eu = 80 m; lower panels, "CASE 2").The shaded regions in (a-d) indicate the shape of the depth-variable diapycnal eddy diffusivity κ (pink), the shape of the irradiance profile (yellow) and the remineralization curve (green) used in the model.The modeled σ t -NO 3 relationship is shown for (e) CASE 1 and (f) CASE 2 and can be compared to fitted profiles from the APEX-ISUS floats from the late winter of 2010 (g, h).

Figure 9 .
Figure 9.The model-derived a over a range of z ML − z eu with a growth rate (µ o ) of 0.2 d −1 (black curve).The location of the curve depends on the µ o selected: µ = 2 d −1 shifts the curve to the right (gray line to the right), and µ = 0.02 d −1 shifts it to the left (gray line to the left).The thin black lines show the model sensitivity to varying upper layer κ.The stars indicate the value of z ML − z eu and a from CASES 1 and 2.

Figure 10 .
Figure10.The a from the APEX-ISUS NO 3 -σ t fits at (a) Hawaii and (b) Bermuda vs. z ML − z eu for each profile (black circles).The method for obtaining approximate z ML and z eu at each profile is described in the Appendix.These observed a are overlaid on the model prediction from Fig.9.
At the time series stations BATS and HOT, the CTD-based profiles of Ch a fluorescence are calibrated to in situ samples of pigment concentration.MB established a statistical relationship between surface-normalized photosynthetically available radiation (PAR) and the depth-integrated Ch a content ( Ch a ) in the water column for Case 1 waters.The relationship between z eu MB and Ch a was approximated with the curves z eu MB = −568.2Chl a −0.746 for z eu MB > −102 m −200.0Chl a −0.293 for z eu MB < −102 m