A comparison of the variability of biological nutrients against depth and potential density

The main biogeochemical nutrient distributions, along with ambient ocean temperature and the light field, control ocean biological productivity. Observations of nutrients are much sparser than physical observations of temperature and salinity, yet it is critical to validate biogeochemical models against these sparse observations if we are to successfully model biological variability and trends. Here we use data from the Bermuda Atlantic Time-series Study and the World Ocean Database 2005 to demonstrate quantitatively that over the entire globe a significant fraction of the temporal variability of phosphate, silicate and nitrate within the oceans is correlated with water density. The temporal variability of these nutrients as a function of depth is almost always greater than as a function of potential density, with he largest reductions in variability found within the main pycnocline. The greater nutrient variability as a function of depth occurs when dynamical processes vertically displace nutrient and density fields together on shorter timescales than biological adjustments. These results show that dynamical processes can have a significant impact on the instantaneous nutrient distributions. These processes must therefore be considered when modeling biogeochemical systems, when comparing such models with observations, or when assimilating data into such models.


Introduction
The distribution of biological nutrients within the world's oceans is one of the significant determining factors in the distribution of oceanic life (Falkowski et al., 1998).For instance, high nutrient coastal regions are the source of most of the world's fish (Watson et al., 2004), while the nutrient poor Correspondence to: J.While (james.while@metoffice.gov.uk)subtropical gyres are relatively devoid of life (Sharp et al., 1980).Consequently, our understanding of biogeochemical processes is intimately linked to our knowledge of the behavior of dissolved nutrients.In particular, modeling efforts (see, for example, Palmer and Totterdell, 2001;Le Quéré et al., 2005) of global ocean biogeochemistry need to accurately describe nutrient fields in order to simulate ocean biogeochemistry.It has long been known that nutrients may exhibit relationships to both ocean temperature (Iselin, 1939;Switzer et al., 2003;Kamykowski, 2008;Steinhoff et al., 2010) and water density (Kamykowski and Zentra, 1985;McGillicuddy Jr. et al., 1999;Archer et al., 1996), with temperature relationships often used to infer nutrient distributions from temperature observations.The form of these relationships locally depends upon surface nutrient availability where the water is formed, and on remineralization and mixing below the surface (Sarimento and Gruber, 2006).Thus the robustness of nutrient-density relationships is determined by many non-local processes.Near to the surface biological consumption is also important, with Gose et al. (2000) showing that surface relationships between nutrients and temperature can be improved by also considering chlorophyll.These previous investigations mostly focus on understanding spatial distributions of nutrients.However, in any location nutrients are also subject to ocean dynamical processes, irrespective of the origin of the water masses.This dynamic temporal variability will always modulate the nutrient variability and determine a component of the observed variations.It is this temporal variability that is the main subject of this study.
Our main interest is in comparing the variability of nutrients about nutrient-potential density (N (ρ)) and nutrientdepth (N (z)) relationships.This is because ocean dynamical processes, particularly the propagation of internal waves and Rossby waves, that vertically displace Lagrangian water properties should have limited impact on N (ρ) variability, but would significantly affect the variability of N (z).Other than ocean dynamics N (z) is determined by much the same processes as N (ρ), although as sinking and remineralization rates are only weakly influenced by changes in potential density, remineralization of sinking particles can be expected to be more strongly correlated to depth.Thus, by comparing the relative variabilities around N(ρ) and N(z), we are able to quantitatively ascertain how much of an impact ocean dynamics has on the instantaneous distribution of nutrients.An analogy to this work can be made by looking at salinity and temperature and similar diagnostics to those used here have been presented for looking at salinity-temperature variability in Troccoli and Haines (1999) and Haines et al. (2006).
Section 2 of this paper shows, in agreement with McGillicuddy Jr. et al. (1999), how the nutrients nitrate, phosphate, and silicate vary within a water column using data taken from the Bermuda Atlantic Time-series Study, (BATS; Phillips and Joyce, 2007).From these abundant data we develop our diagnostic approach to show that the temporal variability of N (ρ) is considerably reduced within the main thermocline.Then in Sect. 3 we develop diagnostics for the much sparser data in the World Ocean Database (WOD05; Garcia et al., 2006) to show that N(ρ) has lower variability than N (z) over almost all of the world's oceans.The paper concludes with a section discussing our results and their relevance to biogeochemical modeling and data assimilation.

Nutrient time-series at the BATS site
To illustrate how the variability of nutrients changes down a water column at a single site, we repeat some of the work of McGillicuddy Jr. et al. (1999), but including more recently available data.Specifically we present results obtained from an analysis of the BATS data set.The BATS time-series, obtained from approximately 75 km SE of Bermuda, provides a long data record of nutrients that stretches back to 1988, with a sampling interval of approximately one month.As such, the BATS data set contains a high density of measurements in a highly localized area.Crucially for our purposes, measurements taken within BATS include phosphate, nitrate (combined with nitrite), silicate, and potential density (derived from temperature and salinity).Such a large collection of data allows us, with a high degree of accuracy, to quantitatively assess the variability of both N(ρ) and N(z).
Nitrate + nitrite, phosphate, and silicate are plotted against depth and potential density (referenced to 2000 m) in Fig. 1 rows (a) and (b).A visual inspection clearly shows the mean relationships between nutrients and depth or density; however, a closer inspection reveals that the scatter of the data around the mean relationship appears much smaller on the N(ρ) plots, (b), than on the N(z) plots, (a).This reduction is most obvious where the nutrient concentration increases dramatically with increasing density through the main thermocline.We quantize this reduction in (c), which shows the variance of N (ρ) and N(z) as a function of depth over the top 1500 m of the water column.The variance of N(z) was cal-culated at 20 m intervals using all data within ±10 m; any interval containing no data was excluded from the result.Conversely, the variance against potential density was obtained by finding the mean density within each 20 m interval and reallocating the nutrient data between the intervals using a nearest in density approach.After this process the nutrient variance was calculated for each depth interval, defined now by its mean density.
Outliers within the data, several of which are clearly visible in Fig. 1b, can severely bias the above variance calculation.As such, outliers were identified and removed prior to calculating the variances.As temperature, salinity, and nutrient errors all contribute to the N (ρ) plots, outliers were identified using these data rather than N (z).Outlier identification was performed by calculating the median absolute deviation of the nutrient data in an advancing 0.4 kg/m 3 window that moved down the water column; any data point lying more than 3 median deviations from the median was flagged as an outlier.Once flagged, outliers were rejected from the variance calculations of both N (ρ) and N (z).
It is clear from the plots in Fig. 1c that between 300 m and 1000 m there is a very substantial reduction in the variability of all of the nutrients.This is indicative of the nutrient distribution being strongly tied to the vertical density structure of the water within the main thermocline.Here dynamically induced variability associated with vertical heaving of the water column raises the depth level variance, but not the density level variance.At other depths, including depths below 1500 m which we do not show, the variance of N (ρ) is approximately equal to the variance of N (z).In shallow waters the biological processing of nutrients and mixing within the mixed layer will tend to break any nutrient-density relationships.In deeper water the dominant variability is less likely to be due to vertical heave of the water column on short timescales and more likely to be due to slower processes.
BATS is but a single time-series and its results apply to only one location.However, it is not unreasonable to assume that similar processes are taking place elsewhere leading to similar patterns of variability with depth.In the next section we extend our scope to look at temporal nutrient variability around the globe.

A global perspective
In order to obtain the variability statistics for nutrients on a global scale, we have analyzed nutrient data collected after 1990 available from the WOD05 database.These data give global coverage of the nutrients phosphate, nitrate, and silicate, though there are large gaps in the data set (the data distribution is shown in Fig. 2a).In order to estimate the variability of N (ρ) and N (z) we rejected any data that were flagged -for whatever reason -within the data set.We also rejected any data that didn't consist of concomitant measurements of temperature, salinity, depth, and nutrient.After these initial checks the data were binned into 2 • × 2 • bins.This bin size was chosen to be large enough to contain sufficient data to allow reasonable statistics to be calculated, yet small enough that lateral variations would be small leaving temporal variability as the dominant signal.Later, by using larger bin sizes, we look at how sensitive the results are to allowing in more spatial variability.

Methods
Outliers in each of these bins were again removed using a median based technique; however, as we have far fewer data in a water column, this process needed to be more flexible than for the BATS data.As with BATS, only N(ρ) data were used in the determination of outliers, but once identified outlying points were rejected from both N(ρ) and N(z).To determine whether a data point k was an outlier it was compared to a set of nearby data points within the water column.These comparison points were found by gradually increasing the density window ρ k ±η k until at least 10 data points were found, where ρ k is the density of k, k is the full potential density range (referenced to the depth of k) of all data within the bin, and η = 0.001,0.002,0.003,...,0.1.If there were fewer than 10 data points available in the largest window (±0.1 ), then all data points within this window were used down to a minimum of 4. The selected data were detrended and the median and median deviation calculated.If the nutrient value of data point k exceeded 3 median deviations from the median of the window, then k was rejected as an outlier.This method detected most outliers; however, it was insensitive to erroneous data points with excessively high densities outside the true density range.To avoid this problem we also checked that the potential density of point k was within 3 median deviations of the potential densities of the other data points within the adaptive window.
With the outliers removed the remaining WOD05 data in each bin were used to estimate the ratio R given by where M is the number of data points per bin down to a user chosen maximum depth.As above, N is the amount of a nutrient (either nitrate, silicate, or phosphate) measured at a data point, ρ is the water density, and z is the depth.• bin to estimate R. The plots were generated from World Ocean Database data measured from 1990 onwards.P <1 is the number of bins with R < 1, P >1 is the number of bins with R > 1, while P <0.9 is the number of bins with R < 0.9.A non-linear color scale is used to highlight the structure near the critical value of R = 1.
We use the ratio R because it is insensitive to factors, such as instrument noise and bias, that would be expected to affect the measured nutrient-density and nutrient-depth relationships equally.If R is less than 1 then there is less scatter, hence less variability, in N(ρ) than in N(z) and vice versa if R > 1.The mean absolute deviation of the data was used in the definition of R because it is a more robust estimator of the variability in the presence of outlying data (DeGroot and Schervish, 2002).This was desirable because our method could not be optimized for every bin, and poor estimates of µ z and µ ρ in bins with sparse data could have lead to outlying data points.
The functions µ ρ and µ z , are the 'mean' nutrient profiles with respect to density and depth in each 2 • × 2 • bin.It is evident from the BATS data, Fig. 1, that these functions may have a complex structure within a water column and we do not seek to explain the origin of these functions from the physical and biogeochemical history of the water masses.For a data point k, a local estimate of µ z was calculated from where z k is the depth of k and µ zk is the value µ z takes at k.The parameters m k and c k , specific to k, are the gradient and intercept of a local linear regression about k.The regression was calculated using all data found within the smallest of z k ± η z that, excluding k, contained 10 data points; here z is the full depth range of the data, and η is defined as before.If the largest window (±0.1 z ) contained less then 10 data points, then all data within the window were used down to a minimum of 3 points, otherwise µ zk could not be estimated and k was rejected.Excluding point k from the estimate ensured that µ zk was independent of k.Since measurements in WOD05 tend to be available at standard depths, it is common to have a lot of data at, or near to, a single depth, with no other nearby data.A linear fit is then ill-conditioned and the data mean was used for µ zk .
The calculation of µ ρ proceeded in a similar fashion, but using the potential density referenced to the depth of k.However, at large depths (higher densities) a complication arose because nutrients, especially silicate, can increase very rapidly with potential density, as seen for silicate in Fig. 1a.With such large gradients µ k becomes almost impossible to determine.To avoid this problem the maximum depth of the analysis was set to 2000 m.

Results
The value of R for the 2 • binned WOD05 data down to a maximum depth of 2000 m is shown for phosphate, silicate, and nitrate in Fig. 2b-d.Also given in this figure are the Biogeosciences, 7, 1263-1269, 2010 number of bins P <1 that have R < 1, the number of bins P >1 for which R > 1, and the number of bins P <0.9 in which R < 0.9.While there is a lot of spatial variability in the plots, it is apparent that most bins, by a ratio of 3:1 or more, have R < 1. Translated into relative area, we have R < 1 over 79% of the area for which we have nitrate data, 80% for phosphate data, and 74% for silicate data.More significantly, for all three nutrients R is less than 0.9 (corresponding to the variability of N (ρ) being at least 11% less than the variability of N(z)) more than twice as often as R > 1.In fact, in terms of relative area, R is less than 0.9 over more than half of the area for which we have data (60% for nitrate; 62% for phosphate; 54% for silicate).Thus from the WOD05 data it appears that temporal variations in nutrients are tied to temporal variations in water density over most of the world's oceans.
There is a lot of scatter in the value of R in Fig. 2, both spatially and between the three nutrients.To get an in-depth and detailed analysis of what is happening would require a bin by bin examination of the data, which is beyond the scope of this paper.Broadly speaking, nitrate and phosphate tend to have smaller values of R than silicate.This is a consequence of silicate tending to have much larger, and more difficult to estimate, nutrient-density gradients at depth.Nonetheless, in the data rich North Pacific and North Atlantic R is generally less than 0.9 demonstrating that N(ρ) shows consistently less, though not much less, scatter than N(z).It may be that with more high quality data in these regions we would expect to get results similar to BATS, where R BATS ≈ 0.7.Interestingly, in the southern oceans below 40 • S, R is often fairly low (sometimes less than 0.5) suggesting good N(ρ) relationships; this is despite the somewhat weaker density stratification present at these latitudes.Although data are limited in the southern oceans, a detailed inspection was conducted on a few bins and the results were found to support the idea of lower variability there.
Further to the above we carried out an experiment on WOD05 nitrate data to see the effect that bin size had on our results.To do this we varied the width of the data bins used to calculated R. The values of P >1 /P <1 obtained from this experiment are shown in Fig. 3.As the bin size grows additional data become available per bin allowing for more accurate statistics, but at the same time an increasing amount of lateral variability is introduced into our results.Clearly for N (ρ), the increased lateral variability does not increase the variance to the same extent as for N(z), leading to the reduction of P >1 /P <1 as the bin size increases.Nutrients, like many other ocean tracers, tend to spread laterally along isopycnals (e.g.McDougall, 1984), and therefore lateral gradients along isopycnals are likely to be smaller than along depth surfaces.For a zero bin size lateral gradients are excluded leaving only temporal variations and linearly extrapolating in Fig. 3 gives 0.31 (regression coefficient of 0.95) for a zero bin size, implying that the temporal variance of N (ρ) relative to N (z) is indeed reduced over ∼70% of the world ocean.Two tests were done to check that the results shown in Figs. 2 and 3 are not due to biases in our analysis method.The first test allowed for a broader data window, with z k ± η z allowed to expanded until it contained 40 data points.This test permitted more data and allowed for more accurate regressions, but was also less local and produced larger errors when µ z and µ ρ were strongly non-linear.Applied to WOD05 nitrate data the results from this test were not significantly different to the results obtained above; the new values were P >1 = 1081, P <1 = 3954, and P <0.9 = 2999.Such small changes in the results show that our method is relatively insensitive to the size of the vertical data window and that we can have confidence in the values of R obtained.
In the second test we carried out the following experiment.An estimate of R was obtained by applying our method to WOD05 depth and density data, but with the WOD05 nutrients replaced with uniformly distributed random values from 0 to the largest measured nutrient value.Discounting outlier removal, not done in the test, we expect R > 1 and R < 1 to be equally likely.The test found a small but distinct bias towards R being greater than 1.This bias is explained by the practice of measuring nutrients on, or as close as possible to, standard depth levels, such as the distinct levels seen on Fig. 1a.Consequently we have a very large number of measurements at relatively few depths and very few measurements elsewhere.This enables µ z to be determined very accurately at depths where we have data.Conversely, plotted against density the data are spread out along a trend, and thus it is more difficult to determine µ ρ .This effect is sufficient to bias R. The bias is slight as, excluding bins with less than 1000 points, the mean value of R was 1.01, with a standard deviation of 0.01.That the plots of Fig. 2 show, despite this bias, a very clear signal of R < 1 almost everywhere, strongly suggests that the variability of nutrients on potential density surfaces is indeed reduced.As with the BATS data, we wish to know where in the water column the reduction in variability of N(ρ) takes place.We calculated the absolute deviations about µ z and µ ρ of all data points in the WOD05 data.These deviations were then collected into 100 m depth bins and their averages taken.The results of this test applied to WOD05 phosphate data (results for nitrate are similar, while silicate follows the same pattern, but has less variability reduction in the near surface) can be seen in Fig. 4a.In the top 1000 m there is a clear reduction of up to 10% in the mean deviation of N (ρ), while below this depth the variability is roughly the same between N (z) and N(ρ).There is some indication that the variability of N (ρ) has been slightly increased at depth, but this effect is small and likely due to the bias of the method.The results seen in the figure are consistent with what is seen at BATS, with upper water column variability reduction and then equality of variability at greater depths.This is in agreement with our earlier assertion that an N (ρ) relationship exists where both nutrients and isopycnals are moved up and down by vertical heave of the thermocline.In a global average this appears to happen over a broader range of depths than in the BATS data because of geographic variations in the mean depth of the thermocline.In regions where the thermocline outcrops, or is very shallow, we get reduced variability near the surface, while reduced variability occurs at depth when the thermocline is very deep.
It is interesting to consider the effects that seasonality has on our results.An in-depth examination of the changes in variability across the globe is challenging due to lack of data, especially in the Southern Hemisphere.However, to see if there is a seasonal change in the variability of the data we conducted a similar experiment to the above.In this experiment we concentrated on nitrate data in the Northern hemisphere, where data are relatively abundant, and then split this data into "summer" (all data from April, May, June, July, August, and September) and 'winter' groups (all data from October, November, December, January, February, and March).Statistics for the value of R were then calculated for both of these groups.It was found for the winter data that P >1 =449 and P <1 =1514, while for the summer data P >1 =541 and P <1 =1778.In both cases the fraction of bins for which R > 1 is very similar (∼ 0.3) implying that there is not much seasonal dependence to our results.From this we may infer that the proportional effect of ocean dynamics on nutrient relationships does not vary significantly throughout the year.

Summary and discussion
We have shown, through the analysis of in-situ nutrient (nitrate, phosphate and silicate) measurements, that a significant amount of the temporal variability of nutrient distributions is coupled to the temporal variability of potential density nearly everywhere within the world's ocean.In other words, nutrient distributions show covariability with potential density.
The reduction in variability around N (ρ) relationships compared to N (z) relationships has been demonstrated both at a single location, using the data rich BATS time-series, and, more importantly, in a global sense using data from the WOD05 database.In both data sets the variability of the nutrients as a function of potential density was shown to be significantly reduced.
Our results, which show little seasonal variability, may be explained by the coupling of nutrients to potential density removing variability due to the dynamical effects of wave propagation.In the case of waves the vertical motion of water affects the nutrients and potential density equally.This is not true for other processes such as biological activity and mixing.These processes add variability to both N (ρ) and N (z).Consequently the smaller variability of N (ρ) when compared to N (z), is most likely due to vertical dynamics inflating the N (z) variability.Thus our results indicate that dynamical vertical motion does play a significant role in determining the instantaneous nutrient structure of the ocean.
In most ocean regions we have shown that the largest reduction in the variability of N (ρ) occurs in upper waters, particularly through the main pycnocline, although not necessarily right to the surface.Close to the surface factors such as biological processing break the N (ρ) relationship.At depth nutrient gradients are small so vertical motion, which is also weak in deep waters, has only a limited impact.However, between these two depth domains significant wave motion acts on steep nutrient gradients leading to increased nutrient-depth, but not nutrient-density, temporal variability.
We suggest that N(ρ) relationships, which we have demonstrated here to be of reduced variability when compared to N (z), can be of importance to practitioners of data assimilation in biogeochemical models.Biogeochemical models can be of very different levels of complexity, but they all share a common dependence on the underlying nutrient distributions.In recent studies data assimilation has been used to constrain either the biogeochemical variables themselves (examples include Triantafyllou et al., 2003;Nerger and Gregg, 2007;Hemmings et al., 2008) and/or the physical state of the ocean (Anderson et al., 2000;Eden and Oschlies, 2006).However, in such assimilation schemes nutrients are usually left unconstrained, leading to the breakdown of the nutrient-water mass relationships and a worsening of the modeled nutrient fields.The results presented here indicate that ideally N (ρ) variability should be kept small, even when assimilating data.This should be possible by applying nutrient balancing increments similar to the approach used for salinity in Troccoli and Haines (1999).We propose that this presents the possibility of considerably improving the reproduction of nutrient distributions in biogeochemical models, potentially having a big impact on all areas of model behavior.Experiments using these ideas are ongoing and results will be reported in a separate study.

Fig. 1 .
Fig. 1.BATS nutrient concentration data plotted against (a) depth and (b) potential density anomaly (ρ −1020 kg/m 3 ) referenced to 2000 m.Plots on row (c) show the nutrient-depth variance (solid-line) and nutrient-density variance (dashed line; see text for how this was calculated) within the top 1500m of the water column.

Fig. 2 .
Fig. 2. The value of R, see text, gridded into 2 • bins.The nutrients shown are (b) Nitrate, (c) Phosphate, and (d) Silicate.(a) Shows the number of data points available per 2 •×2 • bin to estimate R. The plots were generated from World Ocean Database data measured from 1990 onwards.P <1 is the number of bins with R < 1, P >1 is the number of bins with R > 1, while P <0.9 is the number of bins with R < 0.9.A non-linear color scale is used to highlight the structure near the critical value of R = 1.

Fig. 4 .
Fig. 4. Mean absolute deviation of WOD05 phosphate down to 2000 m.Solid line: deviation against depth; dashed line: deviation against density.Plot (b) shows the number of data points available to estimate (a).