Extraction of Photosynthesis Parameters from Time Series Measurements of In Situ Production: Bermuda Atlantic Time-Series Study

Computing the vertical structure of primary production in ocean ecosystem models requires information about the vertical distribution of available light, chlorophyll concentration and photosynthesis response parameters. Conversely, given information on vertical structure of chlorophyll and light, we can extract photosynthesis parameters from vertical profiles of primary production measured at sea, as we illustrate here for the Bermuda Atlantic Time-Series Study. The procedure is based on a model of the production profile, which itself depends on the underwater light field. To model the light field, attenuation coefficients were estimated from measured optical profiles using a simple model of exponential decay of photosynthetically-available irradiance with depth, which accounted for 97% of the variance in the measured optical data. With the underwater light climate known, an analytical solution for the production profile was employed to recover photosynthesis parameters by minimizing the residual model error. The recovered parameters were used to model normalized production profiles and normalized watercolumn production. The model explained 95% of the variance in the measured normalized production at depth and 97% of the variance in measured normalized watercolumn production. A shifted Gaussian function was used to model biomass profiles and accounted for 93% of the variance in measured biomass at depth. An analytical solution for watercolumn production with the shifted Gaussian biomass was also tested. With the recovered photosynthesis parameters, maximum instantaneous growth rates were estimated by using a literature value for the carbon-to-chlorophyll ratio in this region of the Atlantic. An exact relationship between the maximum instantaneous growth rate and the daily growth rate in the ocean was derived. It was shown that calculating the growth rate by dividing the production by the carbon-to-chlorophyll ratio is equivalent to calculating it from the ratio of the final to the initial biomass, even when production is time dependent. Finally, the seasonal cycle of the recovered assimilation number at the Bermuda Station was constructed and analysed. The presented approach enables the estimation of photosynthesis parameters and growth rates from measured production profiles with only a few model assumptions, and increases the utility of in situ primary production measurements. The retrieved parameters have direct applications in satellite-based estimates of primary production from ocean-colour data, of which we give an example.


Introduction
Strict experimental quantification of phytoplankton primary production was initiated with the introduction of the 14 C method in 1952 by Steemann Nielsen [1]. Shortly afterwards followed the quantification of primary production by means of mathematical models based on the first principles of phytoplankton physiology [2,3]. Decades later, we find ourselves with global scale maps of phytoplankton production represented by daily carbon uptake at our disposal [4][5][6][7]. Such maps are made possible by combining satellite measurements of chlorophyll concentration with primary production models forced by photosynthetically-available surface radiation [8,9]. Complementary estimates based on satellite phytoplankton carbon and growth rates have also been developed [10]. At their core, chlorophyll-based models are founded on knowledge of the photosynthesis-light relationship, which enables the calculation of primary production.
The standard method for quantifying the photosynthesis-light relationship is by incubating 14 C enriched phytoplankton samples at a set of light intensities and thereby constructing a photosynthetic response curve [11]. The photosynthesis irradiance function is used to represent the results of these experiments and to describe the photosynthetic response to changes in available light [12,13]. This is accomplished by the non-linear fitting of a prescribed photosynthesis irradiance function to the measurements. The information extracted from the procedure establishes the physiological parameters, which determine the optimum choice of photosynthesis irradiance function from amongst the infinite family of possible functions. These are termed photosynthesis parameters [14,15].
Analogously, the parameters can also be estimated from in situ measurements of primary production [16]. These measurements are carried out at sea with 14 C enriched samples in incubation bottles suspended along a line at predefined depths during the incubation period [17,18]. In this case, the measurements correspond to production integrated in time over the period of incubation, and under natural, variable, irradiance conditions. The standard incubation time is the daylength, and the standard total depth interval is the photic zone [14]. The photosynthesis parameters can be extracted from the measured production profile in a similar fashion to that from incubations under controlled light conditions [19]. The former approach, in which the parameters are extracted from experiments under controlled light conditions, was introduced some decades ago, whereas the latter one was introduced more recently [16].
With the aid of a mathematical model for the daily production profile based on a photosynthesis irradiance function, it became possible to describe the photosynthesis-depth relationship using the same parameters as those used to describe the photosynthesis-light relationship [20]. Such descriptions of the production profile enabled the recovery of photosynthesis parameters from measurements of the production profile that were initially designed to estimate watercolumn production [18].
The example of parameter extraction from measured production profiles shows how mathematical modelling enhances the value of observations carried out at sea. Such indirectly-retrieved values of photosynthesis irradiance parameters, when combined with databases of direct measurements [21], could enhance our understanding of the natural variability in these parameters and the factors that influence them, thereby facilitating model implementations [22,23]. Information about photosynthesis parameters has immediate application in remote sensing, since the estimation of primary production from remotely-sensed data also requires the same photosynthesis parameters [24]. Such indirect estimates of photosynthesis parameters are particularly useful for ocean regions where data on photosynthesis parameters are scarce, but where relatively rich data archives of production profile measurements are available.
Important sources of data for extracting photosynthesis parameters are long-term oceanographic time series, such as the Hawaii Ocean Time-Series [25], the Bermuda Atlantic Time-Series Study [26], the Cariaco Ocean Time-Series [27], the Monterey Bay Time-Series [28], and others worldwide [29,30]. In our previous works, the parameter recovery method was tested and proved to be successful on data from the Hawaii Ocean Time-Series [16,20].
In this work, we extend the application of the method to the Bermuda Atlantic Time-Series Study data set. We begin with a brief description of the model and the data set. Then, we extract the optical parameters from the measured irradiance profiles, as these are required to model the underwater light climate and consequently, the production profile. Following this, we fit an analytical solution for daily production to the measured primary production profiles, which allows extraction of the photosynthesis parameters. We then present their distributions and compare the overall precision of our model for the estimation of normalised production measured at particular depths and for the integrated watercolumn production. Next, a shifted Gaussian function is used to model chlorophyll profiles and to calculate watercolumn production using an analytical solution from Kovač et al. [31]. This function is suitable for describing the vertical profiles of biomass with a deep chlorophyll maximum [32,33], which is usually present at this station. Afterwards, we employ the recovered assimilation numbers to estimate the maximum instantaneous phytoplankton growth rates that are used in carbon-based models and discuss the relationship between the maximum instantaneous growth rate and the daily growth rate. Finally, we construct the seasonal cycle for the assimilation number and use it to model the seasonal cycle of watercolumn production. The novelty of the work is two-fold, practical and theoretical. The practical novelty is the application of a previously established method to a new rich data set from a different locality, and the theoretical novelty is the derivation of an explicit expression linking the growth rate to the production profile.

Model of the Production Profile
The instantaneous primary production per cubic metre at depth in the ocean was determined by the amount of available light, the physiological status of the phytoplankton population and, of course, by the phytoplankton concentration itself [2]. Primary production, P (mg C m −3 h −1 ), is expressed as a rate in carbon units [14,17]. As a measure of phytoplankton biomass throughout this paper, hereafter denoted B, we employed the chlorophyll concentration (mg Chl m −3 ), which is the appropriate measure of biomass for chlorophyll-based models. The chlorophyll molecule couples light with the photosynthetic unit; it is ubiquitous in all phytoplankton taxa and is easily measurable, making it a useful operational index of phytoplankton biomass [34]. As a measure of available light, we took the downward irradiance integrated over the visible portion (400-700 nm) of the spectrum, I (W m −2 ), which is the photosynthetically-active radiation (PAR), easily measured at sea by optical profilers, and can also be modelled relatively easily by knowing just the surface light field and the optical properties of the water column [35,36]. Once B and I are known, what remains to model primary production is the relationship between P and I, i.e., the photosynthesis light function.
This relationship between P and I reflects the physiological status of the phytoplankton and is mathematically expressed as where p B I| α B , P B m is the photosynthesis irradiance function [12], with the following photosynthesis parameters: the initial slope, α B (mg C (mg Chl) −1 (W m −2 ) −1 h −1 ), and the assimilation number, P B m (mg C (mg Chl) −1 h −1 ) [15,37,38]. With this relationship at hand we proceeded to modeling the instantaneous production at depth, P(z). For this we needed knowledge of the underwater light field and we employed a simple model: where K is the diffuse attenuation coefficient for downward irradiance and I 0 is the surface irradiance [34]. Surface irradiance naturally changes with time of day, and, as a consequence, so does the instantaneous primary production. However, primary production measurements in the ocean are made over a certain time interval (typically from sunrise till sunset) at a series of predefined depths [8].
As a result, they give an integral measure of daily primary production, labelled P T (z), expressed as the total amount of carbon assimilated by the phytoplankton per cubic metre, at depth z.
To model P T (z) we assumed, following Platt et al. [39], the following time dependence of surface irradiance: where I m 0 is the noon irradiance and D is the daylength (time from sunrise until sunset). Combining (3) with (2) gave us a model for I(z, t). To calculate P T (z), we integrated (1) over D to get where B(z) is time independent [8]. To solve this particular integral, we the mathematical form of p B (I) needed to be specified. Here, we opted for the Platt et al. [40] function: for which the analytic solution for P T (z) is available [20] and reads where I m * = α B I m 0 /P B m and f z I m * e −Kz is a known analytical function [20]: Dividing both sides of (6) by B(z) gives the solution for normalized daily production P B T (z): The relationship of this production model to other models of the production profile, along with an in-depth discussion of their assumptions and limitations, can be found in Kovač et al. [19].
The quantity on the left hand side of this expression can be estimated simply from measurements by dividing the measured daily production, P T (z), by the measured biomass, B(z), giving P B T (z). Please note here that for the measured or known value of a variable or parameter, we used x, whereas for a model variable or parameter we used x. On the right hand side of (8) noon irradiance is shown, which can also be measured ( I m 0 ), and the attenuation coefficient, which can be determined from profiles of measured downward irradiance ( K). Daylength is easily determined, either from direct observation on the day of measurements or from a model. Here we used a model for daylength from Kirk [34].
Equation (8) shows that P B T (z) is a function of depth with two parameters, α B and P B m . With a set of measurements, P B T (z n ), at a sequence of depths, z n = z 1 , z 2 , ..., z N , we can look for the profile of normalized production (8) that best describes the measurements. This is, in essence, a non-linear optimization problem that can be solved; the exact procedure is described in detail in Kovač et al. [16] and Kovač et al. [20]. In short, we treated P B T (z) as a function of photosynthesis parameters, α B and P B m , and found their optimal values that minimized the error between the model and the measurements. The optimal set of photosynthesis parameters is referred to as the recovered photosynthesis parameters. This procedure was implemented for the Bermuda Atlantic Time-Series Study (BATS) data set, as described below.

Bermuda Atlantic Time-Series Study Data Set
The Bermuda Atlantic Time-Series Study (BATS) station is located at 31 • 40 N 64 • 10 W, some 82 km southeast of Bermuda [41] and within the Western subtropical gyre provice of the North Atlantic, as defined by Longhurst et al. [4]. Measurements as part of the BATS (bats.bios.edu) commenced in October 1988, and the Bermuda Bio Optics Project (oceancolor.ucsb.edu/bbop) commenced in January 1992 [26]. In the data archive for this station, we found a total of 285 cruises with available data. The data set was last accessed on 23 November 2017. Of these, measurements of chlorophyll, primary production, optical profiles and surface PAR were made on 151 cruises.
To be more precise, chlorophyll, primary production and underwater optical measurements were available for all of the 151 cruises, whereas surface irradiance measurements were available for only 87 cruises. For the cruises on which surface irradiance measurements were available, we calculated the noon irradiance, I m 0 , from the total received irradiance , I T , as I m 0 = I T π/2D. Surface PAR was converted to Wm −2 using Smith and Morel's procedure [42]. Out of the 151 cruises, the first full set of measurements we used was from 13 February, 1992 and the last was from 20 May, 2012. The sampling depths for the cruises were 0, 20, 40, 60, 80, 100, 120 and 140 m, for both chlorophyll and primary production. Optical profiles were sampled continuously. Surface PAR was measured with a 2 s interval over the course of daytime for the 87 mentioned cruises. The chlorophyll and primary production data are available at http://batsftp.bios.edu/BATS/bottle/, and the optical data are available at ftp://ftp.eri.ucsb.edu/pub/org/oceancolor/BBOP/. Information on the method used to measure chlorophyll concentration can be found at bats.bios.edu/methods/chapter14.pdf, and information on the method for primary production can be found at bats.bios.edu/methods/chapter18.pdf. Detailed information on the profiling spectroradiometers used to measure vertical optical profiles can be found in Allen et al. [43] in the Material and Methods section, along with other references therein. A thorough review of the BATS measurements was provided by Lomas et al. [44].

Determining the Attenuation Coefficient
A prerequisite for the application of the production model is the knowledge of the underwater light field. In our model this is determined by the I m 0 exp(−Kz) term in the solution for the normalized daily production profile (8). As this is a non-spectral model for underwater irradiance, it requires one value of spectrally-integrated K for the entire PAR domain. We estimated it by first constructing the PAR profile from the measured spectrally resolved irradiance profiles with where I(z n , λ s ) is the measured irradiance per unit wavelength at depth z n and wavelength λ s , with ∆λ s representing the spectral range associated with I(z n , λ s ) and S representing the total number of wavelengths. Here, we assumed that measurements were made simultaneously, which is the same as assuming constant surface irradiance during the cast of the underwater optical profiler, and that the profiles were corrected for any variation in the incoming solar radiation during the course of the measurement. With this procedure, we recovered the PAR profile, I(z n ).
With the knowledge of I(z n ), we proceeded to model it as an exponentially decaying function of depth, as stated in (2), starting from the depth of the first measurement, z 1 , up until the last, z N . These depths, z n , are not the same as the measuring depths for production and chlorophyll, and they change from cruise to cruise because of the nature of the underwater light measurements. To extract the optimal value of K, Equation (2) was fitted to the set of points (z n , I(z n )) and the fitted value of K was later used to model the production profile (8). The underlying assumption behind this procedure is the presence of optical homogeneity in the water column; therefore, the effects of vertically-variable biomass and changing the spectral quality of light on K were neglected. This is justified for oligotrophic open oceans with low chlorophyll concentrations, such as the ocean around the BATS station.
This procedure was applied to all of the profiles available, and the average value of K was 0.0476 m −1 with a standard deviation of 0.0106 m −1 . In total, there were 10,052 point measurements of underwater irradiance, I(z n ), integrated spectrally, according to (9), available from all the cruises. Upon comparing the predicted versus the measured irradiances, the coefficient of determination (expressed in percentages) was 97%. This high coefficient of determination justified the application of (2) as the optical model for this station. With the information on K now available for each cruise, we used it in the model of the production profile to extract the photosynthesis parameters.

Photosynthesis Parameters Extraction
With the measured biomass and production profiles complemented by the estimated attenuation coefficients and the measured noon irradiances, we were able to recover the photosynthesis parameters for this station, as described in Section 2. First, we describe some characteristics of the parameter recovery procedure in relation to the available measurements.
The parameter recovery procedure assigned one set of photosynthesis parameters for the phytoplankton to the entire water column [20]. Similar to the assumption of optical homogeneity imposed to extract one value of the attenuation coefficient for the entire water column, we recognized the underlying assumption of physiological homogeneity for the phytoplankton population. This is justified for the mixed layer, but could be questionable for measurements at greater depths, which may still be in the photic zone. However, for this station, the mixed layer depth and the photic zone are of the same order of magnitude in winter [41], implying that one set of parameters can indeed be assigned to the phytoplankton in the photic zone, as it is most likely well mixed. In summer, the mixed layer depth can be as shallow as 10 to 20 m and in this case, only one incubation was undertaken in the mixed layer.
More precisely, the method recovered the best estimate of the assimilation number for the surface layers, and a best estimate of the initial slope for the deep waters. Typically, P B m determined production rates in light-saturated surface waters, and α B in light-limited deep waters. Greater vertical resolution of the production profile would be required to circumvent this limitation and allow vertically-variable photosynthesis parameters to be assigned [16]. Whether it is justified or not to assign one set of parameters can be judged only retrospectively, by calculating model errors. Given this limitation imposed by the vertical sampling intervals, we proceeded only along these lines, and extracted one set of photosynthesis parameters for the entire water column.
Another obstacle to the straightforward application of the method is the lack of surface irradiance measurements on 64 of the cruises. For these cruises, the median noon irradiance of the remaining cruises was used in place of noon irradiance. The drawback of this is that it prevented us from estimating the initial slope for the mentioned 64 cruises. However, this lack of data was compensated by the fact that α B and I m 0 appeared as a product in the solution for the daily normalized production profile (8). Since there can only be one model production profile that minimizes the model error [16], the optimal profile would have one optimal value for the product, α B I m 0 , and one optimal value for P B m for each measured production profile. Therefore, it was possible to estimate the assimilation number for all the cruises, with or without the measured noon irradiance being available. With respect to the estimation of the assimilation number, this can be identified as a self correcting property of the inverse procedure for parameter recovery. As long as the optimization algorithm can find the global minimum, the extracted value of P B m will indeed be the correct one, i.e., the same as would be obtained by using the measured irradiance. However, we found no way to circumvent the lack of surface irradiance measurements to recover the initial slope. This limited the estimation of the initial slope to 87 cruises which had measured surface irradiance.
In addition, out of all the data, we had to exclude some individual measurements as outliers. We based the exclusion on the criterion that the normalized production profile has to be a decreasing function of depth for a model with uniform photosynthesis parameters [16,20]. We also removed the outliers that were outside two standard deviations from the mean normalized production at each depth. These two criteria combined removed 55 measurements out of a total 1208. In addition, upon applying the method we noticed that it failed on 13 cruises that had a "zig-zag" structure in the measured normalized production profiles. For these profiles, the optimization method did not converge to plausible values of photosynthesis parameters, such that we had to exclude them also. This left us with a total of 138 measured profiles, each with 8 points in the vertical, minus the 55 outliers, totalling 1049 measurements out of the entire data set that were used to recover photosynthesis parameters.
Therefore, we estimated 138 values of the assimilation number and 87 values of the initial slope. The histograms of recovered parameters are given in Figure 1. The mean of the assimilation number was 7.3, the median was 6.7 and the standard deviation was 4.2, all expressed in mg C (mg Chl) −1 h −1 . The mean of the initial slope was 0.3, the median was 0.26 and the standard deviation was 0.27, and only a small number of parameter values fell outside these limits. The P B m distribution ranged mainly from 2 to 15 mg C (mg Chl) −1 h −1 , with a peak at 5 mg C (mg Chl) −1 h −1 and few values outside this range. As there were more P B m values estimated, the distribution of P B m is more reliable than the distribution of α B . With time, as more measurements become available it is expected that these distributions will show a more coherent structure. Having recovered the parameter values, we then used them in the analytical solution (8) to calculate the normalized daily production at depth. We calculated the coefficient of determination which scored a high value of 0.95, meaning 95% of variance in the data set was explained by the model. In addition, we compared all the measurements against the analytical solution, and the comparison is presented in Figure 2. The points on the figure represent the normalized production of the measured biomass divided by P B m D, plotted against I m * e −Kz , as described in detail in Kovač et al. [20]. In total there were 1049 points. From the figure, we can see the points fall in line rather well with the analytical function (7). There are some discrepancies, but the errors are distributed evenly on either side of the model curve, implying that there is no systematic bias in the model for the daily production profile for this data set, justifying the usage of one set of photosynthesis parameters for modelling the normalized production profile at this station. The accuracy of the model was comparable, with respect to the coefficient of variation, to that found for the HOT data set-0.98 for P B T (z) for HOT [20] and 0.95 for P B T (z) for BATS (this work). These results are similar to the accuracy of photosynthesis irradiance models applied to production measured under controlled light conditions [13].
where α B and P B m are the estimated parameters for each profile. The r 2 value between the measured normalized production and the modeled normalized production is 0.95. In total, there are 1049 points.

Estimation of Watercolumn Production
With the recovered parameters at hand, watercolumn production was calculated. We did so by employing the analytical solutions for watercolumn production available from the literature. Of particular interest were the canonical solution for watercolumn production [39] and the solution for watercolumn production with the biomass profile described by the shifted Gaussian function [31].
The two solutions differ with respect to the formulation of the biomass profile. The former uses B(z) = B, a vertically constant biomass, whereas the latter uses B 0 is referred to as the background biomass, h is the integral of the biomass under the Gaussian curve, z m is the depth of the biomass maximum and σ is the width of the biomass peak. The term, h/σ √ 2π, is also abbreviated to H and yields the height of the biomass peak.
The canonical solution is valid for a well-mixed water column with uniform biomass, B, and is given as where f (I m * ) is a known function [39], and the remaining assumptions are the same as in Section 2 of this paper. The solution for the shifted Gaussian biomass augments the canonical solution by an additional term which takes into account the excess production due to vertical variation in biomass. This solution is where ∆P Z,T is the additional term, labelled the stratification term, and the exact expression for it is given in Kovač et al. [31]. At the BATS station, the chlorophyll profile has been typically observed to vary with depth [26,41], and therefore, the canonical solution was not valid for the interval below the mixed layer. Viewed temporally, it is expected to hold for periods of the year during which the mixed layer is deep, such as the winter months at the BATS station [41]. However, it was possible to use the canonical solution to calculate the normalized daily watercolumn production, which represents the production per unit chlorophyll biomass, and hence, is independent of any variation in the chlorophyll profile. Therefore, we were able to test the canonical solution (10) in a manner similar to that in which the solution (8) was tested. The results are shown in Figure 3. Again, the points on the figure were obtained as a combination of measurements and estimated parameters, as described in detail in Kovač et al. [20]. Normalized daily watercolumn production was calculated by the trapezoidal rule from the measurements. Here, the coefficient of determination between the measured and modelled normalized daily watercolumn production was high (0.97). A similarly high coefficient of determination of 0.99 was also found at HOT for normalized watercolumn production [20]. Such high coefficients of determination at BATS and HOT demonstrate that this model structure is appropriate for computing watercolumn production.
Following this, we tested the solution for watercolumn production with the biomass profile described by the shifted Gaussian (11). This solution is tailored for variable biomass profiles, such as those that include a deep chlorophyll maximum, which has been observed to be a quasi-permanent feature at the BATS station. It is known that the shifted Gaussian is a suitable function for the description of this class of biomass profiles [32,45]. Here, we fitted the shifted Gaussian function to the measured biomass profiles, and the comparison of the measured versus modelled biomass is presented in Figure 4a. Fitting was accomplished by tuning the parameters of the shifted Gaussian, namely B 0 , h, σ and z m , to best match the measured biomass profile. The conjugate gradient method was used to perform this task [46,47] and the fit of the shifted Gaussian was convergent for each chlorophyll profile. As can be seen from Figure 4a, the match between the model and the measurement was satisfactory. Quantitatively, the coefficient of determination for the entire biomass data set was 0.93. A total of 1049 chlorophyll measurements were compared with the modelled shifted Gaussian profile in the figure. This is comparable to that obtained for a similar comparison for data from HOT, where the coefficient of determination was found to be either 0.87 or 0.98, depending on whether the deepest two depths were included or excluded [31].  The successful description of the biomass profiles by the shifted Gaussian further justified the application of the solution (11) to calculate watercolumn production. The stratification term, ∆P Z,T , depends on h, σ, and z m , and therefore, the accuracy of the whole solution depends on these parameters. For a certain range of parameter values, this solution exhibited convergent behaviour and gave reasonable estimates of watercolumn production. In particular, when the deep chlorophyll maximum was deep and narrow, the solution converged, whereas when it was shallow and wide, it tended to diverge. As a quantitative measure of convergent behaviour, the 3σ rule was applied, following Kovač et al. [31], where the solution was considered to be valid in situations with z m deeper than 3σ. This is indeed the expected case for the open ocean with a deep chlorophyll maximum, and the 3σ condition was satisfied for 116 profiles from the BATS data set. The results are displayed in Figure 4b. In the figure, we can observe some discrepancies, but an overall matching is evident.

Estimation of Growth Rates
Knowing the assimilation number enabled us to estimate the maximum growth rate of phytoplankton. Following Kirchman [48], Marañon [49] and Sathyendranath et al. [50], we defined the maximum instantaneous growth rate, µ m , attainable by a phytoplankton assemblage as where χ is the carbon-to-chlorophyll ratio. To estimate µ m for a given production profile, we needed information on the assimilation number and the carbon-to-chlorophyll ratio. We took the carbon-to-chlorophyll ratio from Marañon [49] for the Western North Atlantic Subtropical Gyre, equal to 146 mg C (mg Chl) −1 . For each cruise, we calculated µ m , and the results are given in Figure 5, expressed as per day growth rate, by multiplying the hourly rate from Equation (12) by the daylength, D.
In fact, using this method, we determined the maximum daily growth rate. The distribution of the maximum daily growth rate was skewed to the right. The mean value was 0.7 d −1 , the median was 0.62 d −1 and the standard deviation was 0.39 d −1 . The minimum value was 0.13 d −1 and the maximum was 1.67 d −1 . For 118 profiles, the maximum daily growth rate was below one per day, and for 33 profiles, it was above one per day.
The maximum instantaneous growth rates calculated here are slightly higher than those reported in Marañon [49]. However, they reported values for different regions of the Atlantic and did not calculate the µ m for our study area, due to a lack of information on P B m . Also, as stated by the author [49], his estimates were lower than those reported by other authors, such as Gasol et al. [51] and Harris [52]. His reasoning in favour of a low growth rate in the oligotrophic ocean was that phytoplankton do not have sufficiently high assimilation numbers to sustain it. As pointed out by Marañon [49], a µ m > 2.5 d −1 with χ in the range of 100-200 mg C (mg Chl) −1 requires P B m > 25 mg C (mg Chl) −1 h −1 , i.e., above the theoretical maximum calculated by Falkowski [53]. Here, did not observe such high assimilation numbers (Figure 1b), but the histogram shows a right-hand tail with values well above 10 to 15 mg C (mg Chl) −1 h −1 , which implies the possibility of high growth rates, at least on occasion. The bulk of our results are consistent with those of Casey et al. [54], who reported growth rates in the Sargasso Sea equal to 0.42 ± 0.17 d −1 , based on nitrogen uptake rates. The reader is referred to Table 2 in Casey et al. [54] for more literature values on the growth rates in the Sargasso Sea, which are all consistent with the ones reported in this work, although obtained by different methods.
The growth rate calculated with (12) is actually an upper bound on the daily growth rate when reported as a daily rate. Expression (12) is valid for constant and saturating light conditions during the entire daylength. It therefore corresponds to the following growth equation: A more plausible approach would be to allow for the temporal dependence of production on variable light conditions and to calculate the daily growth rate accordingly. Therefore, instead of the upper bound, µ m , the daily growth rate, µ, would be used, starting with the following growth equation: with time-dependent production due to varying irradiance, I(z, t) = I m 0 sin (πt/D) exp(−Kz), where both B and I are now recognized as functions of depth and time. With information on the initial B(z, 0) and final B(z, D) biomass, the daily growth rate at depth z can be expressed as However, we did not have information (measurements) on the final biomass to use (15) directly, but we did have information about the daily production at depth P T (z). To use this information, we expressed the final biomass as a function of the initial biomass and daily production. By integrating (14) over the daylength, D, we obtained the following expression for the final biomass [19]: Upon inserting this expression into (15), we obtained the following equation for the daily growth rate: Going one step further and expressing P B T (z) as a function of depth by using solution (8) we obtained: This expression shows that the daily growth rate at any given depth can be calculated from the maximum instantaneous growth rate (12) and the shape of the production profile, set by the f z I m * e −Kz function. According to (18), the depth dependence of the growth rate is identical to the depth dependence of the normalized production profile in the case of a depth-independent carbon-to-chlorophyll ratio, χ. Whereas it is plausible to have a uniform χ in the mixed layer, it is highly unlikely that χ would remain invariant with depth in a stratified ocean. Therefore, to calculate the daily growth rate for a stratified ocean would require having χ as a function of depth χ = χ(z). For a mixed layer, it may be argued that the daily growth rate should be calculated from mixed-layer production, and this has been studied extensively by Platt et al. [55], who gave explicit expressions for the mixed-layer growth rate. For an in-depth account of mixed-layer growth rate, the reader is also referred to [56,57].
Comparing (12) with (18), we recognized in the numerator of both expressions, that P B m was the maximum possible instantaneous production at depth that would occur under saturating light conditions. According to solution (8), this is never realised for the entire daylength (sunrise till sunset), because the function f z (I m * e −Kz ) is always less than unity and goes to unity only at the limit of I m * e −Kz → ∞ [20]. This limit is not reached in the real ocean, which is a direct consequence of the diurnal cycle of irradiance (3). As the irradiance goes from zero at dawn, to a maximum at noon, and back to zero at sunset, normalized production falls below P B m , and normalized daily production at any depth is always less than P B m D. This is why (12) gives an upper bound on the daily growth rate and (18) gives the daily growth rate, which is less than P B m /χ due to f z (I m * e −Kz ) being less than 1. Therefore, we expected the daily growth rate in the ocean to satisfy the condition µ < P B m /χ when expressed on a per hour basis.

Seasonal Cycle
For each Julian day, we calculated the mean value of the parameters that were estimated from all cruise data measured within the interval set by 15 days before and after that Julian day, for all the years for which data were available. In such a manner, we obtained an estimate of the seasonal cycle of average monthly values of parameters. The seasonal cycle of α B was not as pronounced as the seasonal cycle of P B m . This might be because we only had 87 estimates of α B , whereas we had 151 estimates of P B m . The seasonal cycle of average monthly P B m is given in Figure 6 (thin blue curve). The minimum in P B m occurred in February, after which a steep rise led to a maximum in July. After the maximum, a less steep decline followed until the end of the year. To this seasonal cycle, we further fitted a sum of two sine functions superimposed onto a mean ( Figure 6, thick blue curve). The first sine has a one year period and the second has a half-year period, given here explicitly: P B m (j) = 8.8 − 2.0 sin (2πj/365 + 1.0) + 1.0 sin (4πj/365 − 4.4) , (19) in units of mg C (mg Chl) −1 h −1 , where j = 0, 1, 2, ..., 365 is the Julian day. Judging by the figure, the representation of the seasonal cycle by the sum of two sine functions is well suited for the assimilation number. We further observed that the seasonal cycle of normalized water column production was similar in shape to the seasonal cycle of the assimilation number ( Figure 6, red curve). However, the similarity of the seasonal cycle of watercolumn production was not observed (Figure 7, grey curve). Upon comparing the two, a noticeable feature in the seasonal cycle of P Z,T was a rise in production starting in January and lasting during February, March and April, which has been well documented in the literature [29,30]. This was not observed in the seasonal cycle of P B Z,T , which had its minimum in February and maximum in July. To account for this rise in daily production during the period of January to April, we looked at the seasonal cycle of total chlorophyll concentration within the euphotic zone, (B Z ) (Figure 6, orange curves), and we observed a similar elevation in B Z as in P Z,T during the same period. Therefore the increase in P Z,T can be attributed to a rise in biomass and not in P B m . This is in line with the theoretical solution (10), according to which normalized watercolumn production does not depend on B, but does depend on P B m . Following (10) we can infer that the initial slope, daylength, attenuation coefficient and noon irradiance also play roles, although perhaps not as significant as that of the assimilation number. Further, by comparing the seasonal cycle of B Z with P B m , we observed that they were nearly in counterphase. Starting in January, biomass increased, but P B m declined and then in April, the opposite occurred, with P B m now increasing and B Z decreasing. P B m rose from March untill July, and more or less declined from August until March. On the other hand, B Z declined from February until October (with a steady period during June, July and August) and then increased from October until March. In the context of the observed annual cycle of the assimilation number (Figure 6a) we see that, although the assimilation number was low during winter, the biomass still managed to rise. These findings are consistent with Evans and Parslow [58] for the seasonal cycle in the North Atlantic. According to these authors, the deepening of the mixed layer in winter, combined with low incoming solar radiation caused a decline in the average light available for photosynthesis in the mixed layer, and consequently, a reduction in the average production per unit biomass. A reduction in normalized production was indeed observed at the BATS station (Figure 6), along with an increase in the mixed-layer depth during winter [26,41].
Another consequence of winter mixed layer deepening is the effect it has on zooplankton that graze on phytoplankton. When the mixed layer deepens and the phytoplankton are diluted, this causes a reduction of food for the zooplankton and their numbers decline. This precondition, occurring before the spring, means that when larger phytoplankton start to increase, the larger zooplankton are unable to catch up with them immediately, allowing a bloom to form. Indeed, according to Sarmiento and Gruber [59], the spring bloom at Bermuda is caused by larger phytoplankton, which are able to break free of the grazing control exerted by larger zooplankton. In addition to this, as reported from observations by Chisholm [60], picoplankton are present all year long at a more or less constant biomass in the North Atlantic near Bermuda. This implies that the grazers of smaller phytoplankton are able to catch up on a potential increase in small phytoplankton, so the smaller phytoplankton are kept in check all year long. In simpler terms, as the succession of phytoplankton proceeds over time, the phytoplankton with higher assimilation rates take over and dominate the population because their grazers are unable to catch up on their growth, which subsequently causes a rise in the assimilation number of the total population. Similar behaviour was also reported by Casey et al. [61] and Wallhead [62] and has been observed in multiple-size-class models [63].
Additionally, temperature could also have an effect on the seasonal cycle of the assimilation number. The seasonal cycle of sea surface temperature at Bermuda has been observed to be well defined, with a minimum in February and a maximum in August [30]. The timings of the minimum and maximum temperatures are well aligned with the minimum and maximum assimilation numbers. The amplitude of the seasonal cycle of sea surface temperature was observed to be about 8 • C [30]. Such an amplitude could also have an effect on the assimilation number, since the enzymatic reactions underlying carbon fixation are affected by temperature [53].
Delineating these factors in more detail and decoupling the functional mechanisms would require more measurements and more complex models. Nonetheless, even the measurements we had available, although limited in spatial and temporal resolution, proved useful for the recovery of the assimilation number and the estimation of its seasonal cycle. This in itself is valuable, because prior to this study, there was no information on the seasonal cycle of the assimilation number at BATS station. With time, as more measurements become available, such estimates of the seasonal cycle will be more reliable and precise, and it will become possible to extract the seasonal cycle of the initial slope α B as well.
In an early compilation of seasonal variation in assimilation number [24], the highest assimilation number reported for the Western sub-tropical gyre was 5 mg C (mg Chl) −1 h −1 , in spring. A more recent compilation [21] revised the mean value for spring slightly upwards, to 5.92, which is quite close to the value that we found for P B m in early spring (March, in Figure 6). However, instead of this being a maximum, the seasonal cycle in Figure 6 shows the value in March of~6 mg C (mg Chl) −1 h −1 to be a minimum, with a maximum of~12 mg C (mg Chl) −1 h −1 being reached in summer. The differences may indicate that sporadic data from ships fails to capture the full extent of seasonal variations in photosynthesis irradiance parameters, such as those revealed here. These differences have immediate implications for the computation of marine primary production based on remotely-sensed chlorophyll concentration combined with photosynthesis irradiance parameters (e.g., Sathyedranath et al. [24]).

Application to Remote Sensing
The estimated seasonal cycle of the assimilation number can be used directly to calculate primary production based on remotely-sensed chlorophyll fields, as we demonstrate here. Our goal was to model the seasonal cycle of watercolumn production with remotely-sensed chlorophyll concentration and the estimated seasonal cycle of the assimilation number. For this simple demonstration, we used the canonical solution (10) as the model for watercolumn production on each Julian day (ignoring the influence of any deep chlorophyll maximum in the interest of keeping the demonstration simple). For α B and K, we used yearly average values. We obtained the information about the chlorophyll concentration from the Ocean Colour Climate Change Initiative [64] data archives [65], which are publicly available at www.oceancolour.org. We used the monthly chlorophyll fields and estimated the seasonal cycle of remotely-sensed chlorophyll in the same manner as the seasonal cycle of the assimilation number. We calculated D and I m 0 based on the model from Kirk [34], assuming clear-sky conditions. We employed the canonical solution on each day with the seasonal cycle of the assimilation number based on the sine approximation. The results are presented in Figure 7, along with the average seasonal cycle of measured watercolumn production. In the figure, the annual cycle of primary production obtained with seasonally-varying assimiliation number is labelled P B m (t), and the corresponding result based on the yearly average value of of the assimilation number is labelled P B m . We can see there was indeed a difference between the two estimates. On average, the difference was 8%, with the maximum difference being as high as 24%. The difference was most pronounced during summer. It is also evident that estimates of watercolumn production made with P B m were generally lower compared with the ones made with P B m (t). The results were in better agreement with the measured watercolumn production during winter, which was expected due to the strong mixing and resultant vertical homogeneity in the biomass, which was assumed in (10). During summer, the agreement was not as good. This may be due to three factors: firstly, remotely-sensed chlorophyll has a minimum in summer, whereas in situ surface chlorophyll does not; secondly, the contribution to watercolumn production due to the deep chlorophyll maximum could be significant; and thirdly, the seasonal variations in α B , cloud cover and the diffuse attenuation coefficient, which are ignored here, could be significant. Here, we see clear evidence that, in parts of the seasonal cycle, the stratification term, ∆P Z,T (11), is important. Nevertheless, this simple demonstration highlights the importance of assigning the correct values of photosynthesis irradiance parameters in computations of primary production using satellite data-if the computations had relied on previous archives of the parameters, the computed primary production would have been significantly lower than the observations.

Discussion
In this work, we built upon the previously established procedure for photosynthesis parameter recovery [16,20] by expanding it to estimate the growth rate based on the recovered photosynthesis parameters. We differentiated between the maximum instantaneous growth rate and the daily growth rate and demonstrated the connection between the two. The maximum instantaneous growth rate is defined as the ratio of the assimilation number to the carbon-to-chlorophyll ratio [50]. As such, it represents the growth of phytoplankton at its maximum capacity, as would occur under saturating light conditions. Due to the daily cycle of irradiance, which causes irradiance to drop below saturation levels, the instantaneous growth rate is, at times, below the maximum instantaneous growth rate. To account for this and to relate the daily growth rate to the maximum instantaneous growth rate, we started with a growth equation for phytoplankton, taking into account variable light conditions (14). We demonstrated that if growth effects are taken into account when calculating daily production, then a direct link between daily production and the daily growth rate emerges (17). In fact, we went a step further by providing a direct link between the maximum instantaneous growth rate and the daily growth rate (18). This was done on the premise of Equation (15), in which the daily growth rate is calculated from the logarithm of the ratio of the final to the initial biomass, divided by the time period over which it is calculated. However, some authors take a different approach.
For example, Kirchman [48] argued for not calculating the growth rate by using the expression (15), instead opting to calculate the growth rate by dividing production by the standing stock (biomass). When using (15), he employed B(z, D) = B(z, 0) + P T (z)/χ, which is, in essence, the first order approximation to (16). He then went on to argue that employing (15) is correct only for low growth rates and short incubation times, in other words, the conditions for which the first order approximation holds. Here, we used a solution (16) that takes into account the positive feedback between production and biomass and is valid for any growth rate. When using (16) in (15), the daily growth rate turns out to be calculable as the ratio of production and biomass (17), which is, in essence, what Kirchman [48] advocates as the proper way of calculating the daily growth rate.
Therefore, we have demonstrated that no difference exists between growth rates calculated simply by dividing the production by the standing stock and by calculating the growth rate with (15), when the effect of growth on production is accounted for. Kirchman's (2002) argument for dismissing (15) as the incorrect approach stemmed from using the incorrect expression for the final biomass when calculating the daily growth rate. This is, indeed, not an argument against (15); rather, it is an argument against using B(z, D) = B(z, 0) + P T (z)/χ in (15). If it is acknowledged that phytoplankton increase due to production, then expression (16) for the final biomass should be used in (15).
Another detail to consider is the way in which normalized production is calculated. When using field data, normalized daily production, P B T (z), is calculated by normalizing P T (z) to the initial biomass, thereby obtaining Obviously this problem does not arise in models, because normalized production is calculated with the photosynthesis irradiance function (5). To be more exact, we can use the solution for B(z, D) to express the daily production, P T (z), in the case of the growing biomass given by (14), as in [19]: Combining the previous two expressions gives us the measured normalized daily production on the left hand side and the model normalized daily production on the right hand side: Therefore we see that what we measure does not, in fact, correspond to the true normalized daily production of the model if we allow for growth in biomass during the course of the period of incubation.
Consequently, if P B T (z) deviates significantly from the model, P B T (z), a potential mismatch in the calculated daily growth rates might arise. In the case of the BATS dataset, the model-estimated production with the recovered parameters (α B and P B m ) does not stray far from the measurements (Figure 2), and this implies the calculated growth rates will be approximately the same, whether calculated from the data or from the model-derived expression (18). However, this does come with a penalty, as the assimilation number extracted from (8) might overestimate the true assimilation number estimated from a growth model in which computed production accounted for the growing biomass. The solution for production in this case is given by (21), and we see that solution (6) is a first-order approximation to it (see details in Kovač et al. [19]). Further theoretical work is needed to expand on these considerations.
Another point worthy of discussion is the lack of sufficient data on the carbon-to-chlorophyll ratio which is required in expression (18) for each depth at which one wants to calculate µ(z) from P B T (z), which obliged us to use a fixed value obtained from the literature to calculate the upper bound on the daily growth rate ( Figure 5). An alternative approach would have been to use a photoacclimation model. In the model of Geider et al. [66], with the analytical solution presented by Jackson et al. [67], χ is an explicit function of irradiance scaled to I k : where θ = 1/χ, with θ m being the maximum chlorophyll-to-carbon ratio, and I * = I/I k . In the case of in situ incubations, phytoplankton are maintained at a fixed depth and if we assume adaptation to average daily irradiance at that depth, then we have We can now solve (23) for χ and use it directly in (18), but this only switches the problem to the lack of information on the maximum chlorophyll-to-carbon ratio θ m . However, this is a problem of practice, not of principle, and (18) could be used to calculate the daily growth rate at any depth using the retrieved I k , if θ m were known. Thanks to the parameter recovery procedure, we are now one step closer to estimating the carbon-to-chlorophyll ratio and the growth rates straight from production data measured at sea. The only information still beyond our grasp is θ m , which is required in expression (23).
A related issue is the need for to measure surface irradiance concurrently with in situ primary production measurements. These were unavailable on some cruises and consequently, we were only able to estimate the initial slope for a subset of the cruises, such that there was insufficient information to calculate a reliable parameter distribution and the seasonal cycle for both the initial slope and the photoadaptation parameter. Perhaps with time, as more complete data sets become available, the initial slope will be estimated on more cruises and the seasonal cycle will be better constrained. We might argue that it would be cheaper just to measure surface irradiance in addition to biomass and production profiles than it would be to conduct more laboratory experiments on the photosynthesis light relationship on top of measuring the production profiles to allow estimation of the initial slope. Therefore, from an economic standpoint, measuring surface irradiance routinely has definite advantages.
The surface irradiance also relates to the growth rate. When deriving (17) and (18), we assumed a sinusoidal cycle of daily irradiance. This cycle is an idealisation suitable for a clear sky, but under an overcast sky, the daily cycle changes and the total irradiance reaching the ocean surface is reduced. As a consequence, production is lower in comparison with clear sky irradiance and consequently, the daily growth rate is also lower. Paradoxically, not having an estimate of surface irradiance does not prevent us from estimating the daily growth rate directly, as we can estimate P B m and I m * , and only these quantities appear in (18). Also, we can take measured normalized production directly in (17) and get the growth rate if we know χ. Again, calculating the growth rate precisely is hampered by the lack of information on the carbon-to-chlorophyll ratio.
We should emphasise here the difference between the type of parametrization of phytoplankton photosynthesis and the different types of parametrizations used in atmospheric and physical oceanographic models, which are, in essence, models based on physical laws. In these physical models, parametrizations arise in response to the inability of the models to resolve sub-grid processes [68]. Parametrizations of this sort become superfluous with increased model resolution, as the finer scale models can resolve the previously parametrized physical processes. This stands in sharp contrast to the parametrization of phytoplankton photosynthesis. Here, the parametrization of photosynthesis is itself the fundamental law and is necessary, irrespective of how fine or coarse the model grid spacing is. In essence, the parametrization of photosynthesis cannot be dispensed with on the grounds of increased resolution. This is true not only for the model used here, but for the this whole class of primary production models. In technical terms, the model used is a canonical model, first introduced in the literature by Platt et al. [39]. Strictly speaking, a model is canonical if there is piecewise continuous transformation of variables that transforms any model into the canonical one [69]. This canonical model has the basic structure required to study the entire family of primary production models. For a more detailed description and discussion of the relationship between different models, the reader is referred to Kovač et al. [19] and Kovač et al. [31].
As a final note, we compare our results with those reported earlier. The assimilation number, P B m , reported for the Sargasso Sea by Tin et al. [9] (their Table 2) ranged from 0.88 ± 0.60 up to 7.46 ± 0.33 mg C (mg Chl) −1 h −1 . The same authors [9] also reported that the initial slope, α B , ranged from 0.017 to 0.806 ± 0.46 mg C (mg Chl) −1 (W m −2 ) −1 h −1 . These ranges coincide well with the distribution given in Figure 1. The values of P B m = 7.41 ± 2.09 mg C (mg Chl) −1 h −1 reported by Forget et al. [70] for the Sargasso Sea are also in line with the values shown in Figure 1b. However, the values reported by the same authors for α B = 0.097 ± 0.024 mg C (mg Chl) −1 (W m −2 ) −1 h −1 fall on the lower end of values given in Figure 1a. The numerous values for α B and P B m reported by Platt et al. [71] for the Northwest Atlantic obtained by photosynthesis light experiments all match well with the ones reported in this work. Platt et al. [71] reported that P B m went from about 1 to 12 mg C (mg Chl) −1 h −1 and α B from roughly 0.01 to 0.2 mg C (mg Chl) −1 (W m −2 ) −1 h −1 (the reader is referred to Figure 3a in Platt et al. [71]). We conclude that the values recovered here from in situ production measurements are in a reasonable range for the open ocean and are consistent with previous, direct estimates based on photosynthesis irradiance measurements at sea.
These results have direct implications for the remote sensing of primary production from satellite data.
In comparison, Sathyedranath et al. [24] assigned α B as 0.06 to 0.11 mg C (mg Chl) −1 (W m −2 ) −1 h −1 and P B m as 1.7 to 5.0 mg C (mg Chl) −1 h −1 for this region of the Atlantic. According to the results presented here (Figure 1), these P B m values fall in the lower end of the parameter distributions. Furthermore, rather than applying fixed values for the assimilation number for each season, our results ( Figure 7) allowed a smoothly-varying function to be applied to account for seasonal variation in the parameter when computing primary production using satellite data. These considerations further emphasise the importance of being able to evaluate the seasonal cycle in α B as well.

Conclusions
The presented work dealt with the estimation of photosynthesis parameters from primary production data collected at BATS over the past 30 years. This region of the Atlantic has been well studied experimentally, with a rich archive of in situ primary production measurements [41], but data on photosynthesis parameters are scarce. This is a consequence of the strategies adopted by research institutions [72], which have opted either to measure primary production in situ, or under controlled light conditions, seldom performing both experiments simultaneously. As a result, some time series stations, such as BATS, have a rich data set of in situ primary production measurements, but only a limited data set of photosynthesis parameters.
The parameter recovery procedure applied here established an archive of photosynthesis parameters for BATS ( Figure 1). The recovered parameters were firstly used to model primary production and to test the adequacy of this model in accurately predicting primary production. The overall accuracy in terms of the coefficient of determination was 0.95 for production at depth and 0.97 for watercolumn production. Also, the chlorophyll-depth relation was described with the shifted Gaussian function, which accounted for 0.93 of variance in the measured chlorophyll profiles. Together, such high coefficients of determination demonstrate that this model can accurately calculate primary production at BATS.
The recovered values of the assimilation number were further used to estimate the maximum potential growth rate ( Figure 5), which was previously unknown for this region due to a lack of information on the assimilation number [49]. A mathematical relationship between the maximum instantaneous growth rate and the daily growth rate was also derived (18). A lack of sufficient data on the carbon-to-chlorophyll ratio prevented a more in-depth calculation of the daily growth rates based on the measured production at each depth. Linking the carbon-to-chlorophyll ratio to the growth rate and how it relates to primary production measurements is a potential course for future research.
The seasonal cycle of the assimilation number was estimated, and it was observed to be nearly in counterphase with the chlorophyll seasonal cycle. The highest values of the assimilation number were observed during June, July and August, while the lowest were observed during January, February and March. It was also demonstrated that the seasonal cycle of the assimilation number can be well described by a sum of two sine functions. The utility of this description is that it can easily be used in remote sensing applications to model primary production at this station, as was demonstrated in Figure 7.
Author Contributions: All authors conceived and designed the study. M.L. assisted in acquiring the data. Ž.K. implemented the model and wrote the draft. Ž.K., T.P., S.S. and M.L. wrote the final manuscript. All authors have read and approved the submitted manuscript.