Modeling Water Clarity and Light Quality in Oceans

Phytoplankton is a primary producer of organic compounds, and it forms the base of the food chain in ocean waters. The concentration of phytoplankton in the water column controls water clarity and the amount and quality of light that penetrates through it. The availability of adequate light intensity is a major factor in the health of algae and phytoplankton. There is a strong negative coupling between light intensity and phytoplankton concentration (e.g., through self-shading by the cells), which reduces available light and in return affects the growth rate of the cells. Proper modeling of this coupling is essential to understand primary productivity in the oceans. This paper provides the methodology to model light intensity in the water column, which can be included in relevant water quality models. The methodology implements relationships from bio-optical models, which use phytoplankton chlorophyll a (chl-a) concentration as a surrogate for light attenuation, including absorption and scattering by other attenuators. The presented mathematical methodology estimates the reduction in light intensity due to absorption by pure seawater, chl-a pigment, non-algae particles (NAPs) and colored dissolved organic matter (CDOM), as well as backscattering by pure seawater, phytoplankton particles and NAPs. The methods presented facilitate the prediction of the effects of various environmental and management scenarios (e.g., global warming, altered precipitation patterns, greenhouse gases) on the wellbeing of phytoplankton communities in the oceans as temperature-driven chl-a changes take place.


Introduction
The major factors affecting phytoplankton metabolism are nutrient availability, light and water temperature [1].The maximum depth at which light intensity is adequate to maintain the plant is referred to as the depth of the photic (or euphotic) zone, which is referred to here as the "photic depth" [2].Light intensity at this depth reaches 1% of its surface daylight value, which is enough for photosynthesis to sustain phytoplankton growth and reproduction.This depth can change as the incident solar irradiance changes with time during the day and throughout the year.The photic depth can also change in space as the concentrations of the various attenuators above it change.
The importance of light to phytoplankton has driven many researchers to develop quantitative methods to calculate light intensity in the water column.The complexity of this task arises from the variability not only in the spectral light intensity, but also in the water column due to turbidity.Bio-optical models are used to develop mathematical methods to calculate irradiance through the water column [3].Optical models have to be calibrated based on field measurements of both water turbidity from phytoplankton, colored dissolved organic matter (CDOM) (also called yellow substances or Gelbstoff) and non-algae (non-pigmented) particles (NAPs), as well as associated spectral irradiance [4][5][6][7][8][9].Abdelrhman [10] applied these methods to estuarine systems where both phytoplankton and total suspended solids (TSS) contribute to turbidity.In oceans, the major contributor to turbidity is phytoplankton [11,12].The theoretical basis presented in [12] was further simplified here to accommodate numerical modelers' needs.Prieur and Sathyndranath [13] presented the optical classification of coastal and oceanic waters based on the specific spectral absorption curves of phytoplankton pigments, organic matter and other attenuators of light in the water column.A wide range of absorption and backscattering spectra for oceans was presented by the International Ocean Color Coordination Group (IOCCG) [3].The IOCCG generated a wide range of synthetic data for use in testing remote sensing algorithms.This range was based on mathematical relationships that used measured phytoplankton concentration in ocean waters (Case 1 waters [14]) as a reference for absorption and backscattering by all other attenuators in the water column including NAPs and CDOM.These relationships are used here to calculate irradiance throughout the water column in the oceans.
The focus of this work is determining the available irradiance profile in the water column through the photic depth.The concentration of suspended solids and phytoplankton in the water column controls the amount of light that travels through it.The main objective of this work is to present a mathematical model, which can be included in numerical models, to resolve the coupling between irradiance and phytoplankton in the oceans.This objective is met by focusing on using available information from bio-optical methods, rather than developing them.Inherent optical properties of water (i.e., chl-a concentrations) [3] are used to meet the modeling objective.

Methods
The total concentration of chl-a is obtained from available field measurements and used to develop the mathematical methodology to estimate the irradiance throughout the water column at any location throughout the whole year.Figure 1 presents a definition sketch of the vertical structure for the mathematical model for irradiance, and Table 1 presents the definitions of all abbreviations and symbols.Attenuation of the incident irradiance through the water column is a function of absorption and scattering by the pure seawater and the dissolved and particulate materials therein.
contributor to turbidity is phytoplankton [11,12].The theoretical basis presented in [12] was further simplified here to accommodate numerical modelers' needs.Prieur and Sathyndranath [13] presented the optical classification of coastal and oceanic waters based on the specific spectral absorption curves of phytoplankton pigments, organic matter and other attenuators of light in the water column.A wide range of absorption and backscattering spectra for oceans was presented by the International Ocean Color Coordination Group (IOCCG) [3].The IOCCG generated a wide range of synthetic data for use in testing remote sensing algorithms.This range was based on mathematical relationships that used measured phytoplankton concentration in ocean waters (Case 1 waters [14]) as a reference for absorption and backscattering by all other attenuators in the water column including NAPs and CDOM.These relationships are used here to calculate irradiance throughout the water column in the oceans.
The focus of this work is determining the available irradiance profile in the water column through the photic depth.The concentration of suspended solids and phytoplankton in the water column controls the amount of light that travels through it.The main objective of this work is to present a mathematical model, which can be included in numerical models, to resolve the coupling between irradiance and phytoplankton in the oceans.This objective is met by focusing on using available information from bio-optical methods, rather than developing them.Inherent optical properties of water (i.e., chl-a concentrations) [3] are used to meet the modeling objective.

Methods
The total concentration of chl-a is obtained from available field measurements and used to develop the mathematical methodology to estimate the irradiance throughout the water column at any location throughout the whole year.Figure 1 presents a definition sketch of the vertical structure for the mathematical model for irradiance, and Table 1 presents the definitions of all abbreviations and symbols.Attenuation of the incident irradiance through the water column is a function of absorption and scattering by the pure seawater and the dissolved and particulate materials therein.
Exponentiation base (e = 2.718281828459) Calibration coefficients dimensionless RR i (λ j ) Reduction ratio of wavelength λ j at the i-th layer dimensionless S NAPs concentration mg•L −1 SS Spectral slope nm −1 λ Wavelength nm

Subscripts i
Counter for the water layer, j Counter for the wavelength λ t Time x Eastward spatial location (see Figure 1) y Northward spatial location (see Figure 1) z Vertical spatial location below the water surface (see Figure 1)

+
Normalized value (see Table 2 and Figure 2) a To avoid confusion in the subscripts presented in various equations, b is used for backscattering instead of the commonly used b b , and K is used instead of K d for the downwelling extinction coefficient.
There are seven major contributors to the loss of light intensity through the water column:  (7) backscattering by pure seawater [3,8,13,15].While volume scattering exists in all directions, only backscattering is usually considered as a loss in the downwelling irradiance [2].Although some methods combine some of these basic contributors, the following methodology provides calculations of these seven types of losses within the range of the photosynthetically-available radiation (PAR) (400 nm-700 nm).
The following equations define the mathematical values for the absorption and backscattering coefficients according to the basic relationships presented in [3] for Case 1 waters.These relationships use phytoplankton chl-a concentration (µg•L −1 = mg•m −3 ) as a reference for the absorption relationships of phytoplankton pigment, NAPs and CDOM, as well as backscattering from phytoplankton particles and NAPs.Light absorption by phytoplankton, NAPs and CDOM uses reference absorption values of phytoplankton pigment at the wavelength λ = 440 nm, while light backscattering uses reference phytoplankton backscattering values at λ = 550 nm.Absorption coefficients for the whole visible range of the spectrum are calculated using the normalized spectral absorption values provided in the literature.Figure 2 and Table 2 present spectral distributions related to absorption by pure seawater, phytoplankton pigment, NAPs and CDOM; in addition to spectral distributions related to backscattering from phytoplankton, NAPs and seawater, as well as the derived spectral distribution function for the incident light [10].
The following equations define the mathematical values for the absorption and backscattering coefficients according to the basic relationships presented in [3] for Case 1 waters.These relationships use phytoplankton chl-a concentration (μg•L −1 = mg•m −3 ) as a reference for the absorption relationships of phytoplankton pigment, NAPs and CDOM, as well as backscattering from phytoplankton particles and NAPs.Light absorption by phytoplankton, NAPs and CDOM uses reference absorption values of phytoplankton pigment at the wavelength λ = 440 nm, while light backscattering uses reference phytoplankton backscattering values at λ = 550 nm.Absorption coefficients for the whole visible range of the spectrum are calculated using the normalized spectral absorption values provided in the literature.Figure 2 and Table 2 present spectral distributions related to absorption by pure seawater, phytoplankton pigment, NAPs and CDOM; in addition to spectral distributions related to backscattering from phytoplankton, NAPs and seawater, as well as the derived spectral distribution function for the incident light [10].2 for the definitions of normalized absorption and backscattering coefficients.Table 2. Spectral distributions of absorption, a, backscattering, b, and the shape function, f, for the PAR (see the footnotes for details).2 for the definitions of normalized absorption and backscattering coefficients.

Overall Light Attenuation
The mathematical structure of the model to calculate irradiance in oceans is similar to that presented by [10] for estuaries.The only difference is the photic depth, which defines the lower bound for irradiance in the oceans rather than the bed in shallower systems.The intensity of light at any depth can be represented by the Beer-Lambert law.The following equations summarize this model.
Beer-Lambert law: Spectral irradiance at the bottom of the first (surface) layer: with the incident irradiance: and the spectral extinction coefficient: Irradiance at the bottom of a general layer, i: Overall irradiance for numerical integration: where: where E 0 is the incident irradiance just beneath the water surface is the extinction coefficient of the downwelling spectral irradiance (m −1 ) and f (λ) is the distribution function of the incident light between the various wavelengths within the PAR [10].The subscripts w, c, s, g and p refer to water, chl-a, NAPs, CDOM (Gelbstoff) and phytoplankton, respectively; and the superscript i indicates the layer number.RR i (λ j ) is the reduction ratio (dimensionless) of the incident light within λ j through layer i .The subscript j refers to the discrete values of the normalized absorption coefficients at λ j values representing PAR in 5-nm increments (j = 1-61; Table 2).According to Simpson's rule, only half of the first and last values can be used in each summation (i.e., at j = 1 and j = 61).Introducing the 5-nm increment in Equation ( 6) preserves the total irradiance within the PAR.The numerical integration procedure is executed for each layer within the water column.The superscripts and subscripts are sometimes dropped for convenience.The same consistent notation is used in the following descriptions of the various attenuators at the same discrete wavelengths considering λ and λ j to be synonymous.

Light Absorption by Pure Seawater, a w
Prieur and Sathyendranath [13] presented the absorption coefficient values, a w (λ j ), for pure seawater at discrete wavelengths, λ j (Table 2, Figure 2).

Light Absorption by Algal Pigment, a c
Light absorption by phytoplankton pigment is given by the following equations [16,17]: a c (440) = 0.05 C x,y,z,t 0.626 (9) where a c (λ) is the absorption coefficient by phytoplankton pigment (m −1 ) at any wavelength λ, a c (440) is the absorption coefficient by phytoplankton pigment (m −1 ) at wavelength 440 nm, a c + (λ) is the normalized spectral absorption value at wavelength λ with respect to absorption at λ = 440 nm [13] and  2) provide a consistent parameterization for the modeling methodology presented here.As expected, at the normalization wavelength, λ = 440 nm, a c + (440) = 1 (Table 2).Gallegos [18] indicated that the values of Prieur and Sathyendranath [13] are adequate for use in his bio-optical methods.The calculated value of phytoplankton absorption, a c (440), is used in the following calculations of the other absorption and backscattering coefficients.

Light Absorption by NAPs, a s
The following equations are from IOCCG [3]: a s (440) = P 1 a c (440) where a s (λ) is the absorption coefficient by NAPs (m −1 ) at any wavelength λ, a s (440) is the absorption coefficient by NAPs (m −1 ) at wavelength 440 nm and R 1 is a random value between 0.0 and 1.0.The randomness in R 1 controls the random values of P 1 makes the relationship between a s (440) and a c (440) not fixed and avoids extremely large a s (440) when a c (440) is small.The range of the random variable P 1 is 0.1-0.6, and its distribution is presented in [3].SS s is the spectral slope for NAPs (randomly valued between 0.007 and 0.015 nm −1 [3]).Recent studies indicate that the NAPs vs. the λ absorption curve has an exponential decay shape [19,20] similar to CDOM.However, understanding of the NAPs behavior is still very limited, and more detailed studies are recommended [20].Until future values become available, this work assumes that the spectral slope for NAPs is in the middle of the above range (i.e., SS s = 0.011 nm −1 ).To eliminate the randomness, R 1 is calibrated as presented in the Calibration and Validation Section.Values of the spectral distribution a + s (λ) = exp(−SS s (λ − 440)) are presented in Table 2 and Figure 2.

Light Absorption by CDOM, a g
The following equations are from IOCCG [3]: where a g (λ) is the absorption coefficient by CDOM (m −1 ) at any wavelength λ, a g (440) is the absorption coefficient by CDOM (m −1 ) at a wavelength of 440 nm, SS g is the spectral slope for CDOM between 0.01 and 0.02 nm −1 and R 2 is a random number between 0.0 and 1.0.The randomness in R 2 controls the random values of P 2 , makes the relationship between a g (440) and a c (440) not fixed and avoids extremely large a g (440) when a c (440) is small.The range of the random P 2 values is 0.3-6.0, and its distribution is presented in [3].To eliminate the randomness, R 2 is calibrated as presented in the Calibration and Validation Section.Values of the spectral distribution a + g (λ) = exp(−SS g (λ − 440)) are presented in Table 2 and Figure 2. b p (550) = P 3 C x,y,z,t 0.57 ( 17) where b p (λ) is the backscattering of phytoplankton at wavelength λ, b p is the backscattering fraction, which depends on the phase function of phytoplankton (assumed 1% based on the Fournier-Forand phase function with respect to scattering angle [21]), P 3 is randomly valued between 0.06 and 0.6 for a given C x,y,z,t and R 3 is a random value between 0.0 and 1.0.The range of the random n 1 values is −0.1-2.0, and its distribution is presented in [3].The randomness in R 3 controls the random values of n 1 , makes the relationship between n 1 and C x,y,z,t not fixed and avoids extremely large n 1 when C x,y,z,t is small.The randomness in P 3 controls the random values of b p (550) and makes the relationship between b p (550) and C x,y,z,t not fixed.The range of the random values and distribution of n  2 and Figure 2.These values will change with depth according to the chl-a concentration (Equation ( 18)).Similarly, the calculations of backscattering for NAP particles include the wavelength-dependent parameters for backscattering for the whole 400-700-nm spectrum, which are based on normalized values referenced to the wavelength λ = 550 nm.The following equations are from IOCCG [3]: b s (550) = P 4 C x,y,z,t 0.766 ( 20) where b s (λ) is the backscattering of NAPs at wavelength λ, b s is the backscattering fraction, which depends on the average particle phase function of phytoplankton (assumed 0.0183 based on the Petzold phase function with respect to scattering angle [22]), P 4 is randomly valued between 0.06 and 0.6 for a given C x,y,z,t and R 4 is a random value between 0.0 and 1.0.The randomness in R 4 controls the random values of n 2 , makes the relationship between n 2 and C x,y,z,t not fixed and avoids extremely large n 2 when C x,y,z,t is small.The randomness in P 4 controls the random values of b s (550) and makes the relationship between b s (550) and C x,y,z,t not fixed.2 and Figure 2.These values will change with depth according to the chl-a concentration (Equation ( 21)).
2.1.7.Light Backscattering by Pure Seawater, b w Buiteveld [23] presented the backscattering coefficient values, b w (λ), for pure seawater at 10-nm increments within the visible range.These values were linearly interpolated at 5-nm increments to fit within the same range of the numerical integration for the other coefficients (Table 2, Figure 2).

Data
Existing irradiance and chl-a data from the North Atlantic were used to calibrate and validate the presented methodology.Data were obtained through personal communication [24].Details of the field methods, protocols and times are published in [25], and they are not repeated here.The data included incident surface irradiance (W•m −2 ) at the water surface for all wavelengths within PAR at 5-nm increments, profiles of the downwelling irradiance (W•m −2 ) at 2 m-depth increments and chl-a concentration profiles at various depth locations during various seasons.The data indicated that chl-a concentrations were always <3.5 µg•L −1 .Data analysis indicated that the ship-based and the glider-based surface PAR agree [25].Accordingly, the glider surface records were used here to represent the incident irradiance for the remainder of the corresponding downwelling profile.
Two sets of consistent data of incident irradiance, irradiance profile and chl-a profile were used to calibrate and validate the presented methodology.The major criteria to choose these data were: cloud cover does not obstruct incident irradiance during the period 10:00 GMT-14:00 GMT; the chl-a profile measurements are collected on the same day of the irradiance data; and the data extend to at least 100 m below the water surface.The collected chl-a concentrations were sparse throughout the monitored top part of the water column and had to be interpolated at the same depth increments as the irradiance profile (i.e., 2 m).Only fourteen chl-a data profiles reached depths ≥100 m.Two of these profiles with contrasting vertical distributions were used here.

The Spectrum Shape Function, f
For modeling purposes, Abdelrhman [10] suggested the use of the spectrum shape function at the top of the atmosphere to redistribute the incident irradiance, E 0 , to its spectral values, E 0 (λ), within the PAR at the water surface (Table 2).The spectrum shape function alleviates the need to identify E 0 (λ) at each location throughout the simulation time.Instead, the overall incident irradiance, E 0 , can be used with f (λ) (Equation (3)).This simplification is tested with the incident irradiance of 339 W•m −2 on 11 June 2013 at 9:30 GMT. Figure 3a compares the shapes of the actual incident irradiance to that of the shape function, f (λ) (Table 2).The difference between the calculated and observed values is ±5% for most of the PAR (λ = 450 nm-650 nm) and <±15% outside that range.The scatter plot (Figure 3b) indicates that the correlation (slope) between calculated and observed shapes is 0.9981 (R 2 = 0.7068) with zero intercept, which indicates that the use of f (λ) is adequate for modeling purposes.The spectrum shape function is used here to define the spectral distribution for the topmost near-surface record of the downwelling irradiance profile.
the PAR at the water surface (Table 2).The spectrum shape function alleviates the need to identify E0(λ) at each location throughout the simulation time.Instead, the overall incident irradiance, E0, can be used with f(λ) (Equation ( 3)).This simplification is tested with the incident irradiance of 339 W•m −2 on 11 June 2013 at 9:30 GMT. Figure 3a compares the shapes of the actual incident irradiance to that of the shape function, f(λ) (Table 2).The difference between the calculated and observed values is ±5% for most of the PAR (λ = 450 nm-650 nm) and <±15% outside that range.The scatter plot (Figure 3b) indicates that the correlation (slope) between calculated and observed shapes is 0.9981 (R 2 = 0.7068) with zero intercept, which indicates that the use of f(λ) is adequate for modeling purposes.The spectrum shape function is used here to define the spectral distribution for the topmost near-surface record of the downwelling irradiance profile.

The Summer Dataset
The summer data used in this study were collected on 11 June 2013 (Figure 4).The concentration of chl-a showed a systematically decreasing trend with depth below the water surface (Figure 4a).This trend was used to calculate chl-a concentrations at 2 m-depth increments, which correspond to the measurements of the downwelling irradiance.The average of four irradiance profiles collected during midday (10:00 GMT-14:00 GMT) was used to represent the observed downwelling irradiance, E, on that day (Figure 4b).The topmost irradiance value was assumed to represent the incident irradiance, E 0 .The distribution of the incident PAR was reconstructed according to Equation (3) (Figure 4c).

The Summer Dataset
The summer data used in this study were collected on 11 June 2013 (Figure 4).The concentration of chl-a showed a systematically decreasing trend with depth below the water surface (Figure 4a).This trend was used to calculate chl-a concentrations at 2 m-depth increments, which correspond to the measurements of the downwelling irradiance.The average of four irradiance profiles collected during midday (10:00 GMT-14:00 GMT) was used to represent the observed downwelling irradiance, E, on that day (Figure 4b).The topmost irradiance value was assumed to represent the incident irradiance, E0.The distribution of the incident PAR was reconstructed according to Equation (3) (Figure 4c).

The Fall Data Set
A subsurface maximum of chl-a concentration appeared in the fall season of 2013 at ~30-40 m below the water surface [25].The data on 3 September 2013 consisted of the chl-a profile (Figure 5a), which demonstrated the subsurface maximum chl-a concentration.Two trends were used to interpolate chl-a concentrations at 2-m increments corresponding to the measured downwelling irradiance.The average of four irradiance profiles collected during midday (10:00 GMT-14:00 GMT) was used to represent the observed downwelling irradiance, E, on that day (Figure 5b).The topmost irradiance value was assumed to represent the incident irradiance, E0.The distribution of the incident

The Fall Data Set
A subsurface maximum of chl-a concentration appeared in the fall season of 2013 at ~30-40 m below the water surface [25].The data on 3 September 2013 consisted of the chl-a profile (Figure 5a), which demonstrated the subsurface maximum chl-a concentration.Two trends were used to interpolate chl-a concentrations at 2-m increments corresponding to the measured downwelling irradiance.The average of four irradiance profiles collected during midday (10:00 GMT-14:00 GMT) was used to represent the observed downwelling irradiance, E, on that day (Figure 5b).The topmost irradiance value was assumed to represent the incident irradiance, E 0 .The distribution of the incident PAR was reconstructed according to Equation (3) (Figure 5c).

Calibration and Validation
The following three important rules of thumb were observed in calibration: (1) any attenuation coefficient for any contributor at any wavelength (i.e., the exponentiation terms in Equations ( 5)-( 7)) cannot exceed 1.0 at any depth; (2) any attenuation coefficient for any contributor except water should be correlated with chl-a concentration (the surrogate for all attenuations); and (3) observed irradiance should not exceed the irradiance in pure seawater.

Calibration with Fall Data
The main purpose of the calibration is to define non-random (fixed) values for R1, R2, R3, R4, P3 and P4.The measured topmost value of the downwelling irradiance represented E0 on 3 September 2013; the recorded light intensity at each depth represented Eℓ; depth increments (2 m) represented the layer thickness, ℓ; and the interpolated chl-a concentration within each layer represented Cx,y,z,t.The calibration started with the mid-range values of the parameter and proceeded to improve the fit with observations of the downwelling irradiance to find the correlation between predicted and observed values.A zero intercept was enforced in all of the correlation plots.For convenience, the p-values were assumed to be the same within their range of 0.06-0.6(range = 0.54), and the R values were also assumed to be the same within their specified range of 0.0-1.0(range = 1.0).Both p and R ranges were divided into 10 increments.Systematic iterations took place by fixing the p-values at each of the incremented values and checking the correlation between predicted and observed irradiances at each of the R incremented values.Figure 6a presents the calibrated irradiance profile with all of the R parameters calibrated to 0.3 and the p parameters calibrated to 0.438.The correlation scatter plot (Figure 6b) indicates that the calculated irradiance agreed with observed values (R 2 = 0.9925).The regression line slope was 0.9994, which is very close to unity for a one-to-one correlation.

Calibration and Validation
The following three important rules of thumb were observed in calibration: (1) any attenuation coefficient for any contributor at any wavelength (i.e., the exponentiation terms in Equations ( 5)-( 7)) cannot exceed 1.0 at any depth; (2) any attenuation coefficient for any contributor except water should be correlated with chl-a concentration (the surrogate for all attenuations); and (3) observed irradiance should not exceed the irradiance in pure seawater.

Calibration with Fall Data
The main purpose of the calibration is to define non-random (fixed) values for R 1 , R 2 , R 3 , R 4 , P 3 and P 4 .The measured topmost value of the downwelling irradiance represented E 0 on 3 September 2013; the recorded light intensity at each depth represented E ; depth increments (2 m) represented the layer thickness, ; and the interpolated chl-a concentration within each layer represented C x,y,z,t .The calibration started with the mid-range values of the parameter and proceeded to improve the fit with observations of the downwelling irradiance to find the correlation between predicted and observed values.A zero intercept was enforced in all of the correlation plots.For convenience, the p-values were assumed to be the same within their range of 0.06-0.6(range = 0.54), and the R values were also assumed to be the same within their specified range of 0.0-1.0(range = 1.0).Both p and R ranges were divided into 10 increments.Systematic iterations took place by fixing the p-values at each of the incremented values and checking the correlation between predicted and observed irradiances at each of the R incremented values.Figure 6a presents the calibrated irradiance profile with all of the R parameters calibrated to 0.3 and the p parameters calibrated to 0.438.The correlation scatter plot (Figure 6b) indicates that the calculated irradiance agreed with observed values (R 2 = 0.9925).The regression line slope was 0.9994, which is very close to unity for a one-to-one correlation.

Validation with Summer Data
Using the same calibration values for the R and p parameters, the summer irradiance profile on 11 June 2013 was checked against observations (Figure 7a).The scatter plot showed weaker linear agreement with the r 2 value of 0.8497 (Figure 7b).The regression line slope with zero intercept was 0.709, which reflects a one to-one correlation between predicted and observed profiles.The reason for the lower correlation was attributed to other factors (e.g., stratification).In addition, the reported irradiance measurements within the top 12 m coincided with the values of pure seawater (Figure 7a), which un-intuitively does not represent effects from other attenuators.Validation with other data is needed (see the discussion).

Results
The following results present a sample to illustrate the availability and quality of light in ocean waters.The availability of light is defined by the photic depth at which 1% of the incident irradiance exists.In addition, the presented methodology can identify the various depths where the spectral values reach 1% of their incident values.As some spectral bands will become extinct before others, the quality of light refers to the amount of energy remaining in each spectral band at any depth.

Validation with Summer Data
Using the same calibration values for the R and p parameters, the summer irradiance profile on 11 June 2013 was checked against observations (Figure 7a).The scatter plot showed weaker linear agreement with the r 2 value of 0.8497 (Figure 7b).The regression line slope with zero intercept was 0.709, which reflects a one to-one correlation between predicted and observed profiles.The reason for the lower correlation was attributed to other factors (e.g., stratification).In addition, the reported irradiance measurements within the top 12 m coincided with the values of pure seawater (Figure 7a), which un-intuitively does not represent effects from other attenuators.Validation with other data is needed (see the discussion).

Validation with Summer Data
Using the same calibration values for the R and p parameters, the summer irradiance profile on 11 June 2013 was checked against observations (Figure 7a).The scatter plot showed weaker linear agreement with the r 2 value of 0.8497 (Figure 7b).The regression line slope with zero intercept was 0.709, which reflects a one to-one correlation between predicted and observed profiles.The reason for the lower correlation was attributed to other factors (e.g., stratification).In addition, the reported irradiance measurements within the top 12 m coincided with the values of pure seawater (Figure 7a), which un-intuitively does not represent effects from other attenuators.Validation with other data is needed (see the discussion).

Results
The following results present a sample to illustrate the availability and quality of light in ocean waters.The availability of light is defined by the photic depth at which 1% of the incident irradiance exists.In addition, the presented methodology can identify the various depths where the spectral values reach 1% of their incident values.As some spectral bands will become extinct before others, the quality of light refers to the amount of energy remaining in each spectral band at any depth.

Results
The following results present a sample to illustrate the availability and quality of light in ocean waters.The availability of light is defined by the photic depth at which 1% of the incident irradiance exists.In addition, the presented methodology can identify the various depths where the spectral values reach 1% of their incident values.As some spectral bands will become extinct before others, the quality of light refers to the amount of energy remaining in each spectral band at any depth.
Figure 8 presents the percent change in the total downwelling irradiance with depth, which indicates a photic depth of 50 m on 3 September 2013.Figure 9 presents the penetration profiles of various wavelengths.The high wavelengths are attenuated heavily within the top 10-20 m, while the shorter wavelengths continue to deeper waters.The penetration depths of the whole PAR on the same day are shown in Figure 10.Wavelengths 485-500 nm reached a maximum depth of 46 m, while those close to 700 nm were almost extinct at an 8-m depth.Figure 11 presents the change in the spectral distribution of the PAR throughout the water column, which defines the light quality on that day.It is worth recalling that the shape and penetration depth of the various wavelengths will change throughout space and time as the chl-a concentration and incident irradiance change.
Figure 8 presents the percent change in the total downwelling irradiance with depth, which indicates a photic depth of 50 m on 3 September 2013.Figure 9 presents the penetration profiles of various wavelengths.The high wavelengths are attenuated heavily within the top 10-20 m, while the shorter wavelengths continue to deeper waters.The penetration depths of the whole PAR on the same day are shown in Figure 10.Wavelengths 485-500 nm reached a maximum depth of 46 m, while those close to 700 nm were almost extinct at an 8-m depth.Figure 11 presents the change in the spectral distribution of the PAR throughout the water column, which defines the light quality on that day.It is worth recalling that the shape and penetration depth of the various wavelengths will change throughout space and time as the chl-a concentration and incident irradiance change.Figure 8 presents the percent change in the total downwelling irradiance with depth, which indicates a photic depth of 50 m on 3 September 2013.Figure 9 presents the penetration profiles of various wavelengths.The high wavelengths are attenuated heavily within the top 10-20 m, while the shorter wavelengths continue to deeper waters.The penetration depths of the whole PAR on the same day are shown in Figure 10.Wavelengths 485-500 nm reached a maximum depth of 46 m, while those close to 700 nm were almost extinct at an 8-m depth.Figure 11 presents the change in the spectral distribution of the PAR throughout the water column, which defines the light quality on that day.It is worth recalling that the shape and penetration depth of the various wavelengths will change throughout space and time as the chl-a concentration and incident irradiance change.

Discussion
A mathematical approach was presented to calculate irradiance in ocean waters (Case 1 waters).The methodology is not site-specific, and it can be applied to any system.Nevertheless, proper calibration of the bio-optical model remains an essential factor for any site-specific application.The methodology presented here utilized relationships from a wide range of absorption and backscattering spectra that were presented by the IOCCG [3].This range was based on mathematical relationships that used phytoplankton concentration in ocean waters as a reference to estimate all absorption and backscattering components in the water column.Employing numerical modeling with bio-optical modeling validates the use of the presented approach for predictions of future scenarios related to changes in environmental and anthropogenic conditions.For example, environmental impacts due to global warming may cause alterations to seasonal cycles in temperature, precipitation and wind patterns.Such alterations can impact ocean circulation and stratification, which directly impacts the transport, distribution and composition of phytoplankton groups and consequently their optical properties.Similarly, anthropogenic effects due to increased aerosols can alter the incident light, which affects the photic depth and the wellbeing of the phytoplankton communities.

Discussion
A mathematical approach was presented to calculate irradiance in ocean waters (Case 1 waters).The methodology is not site-specific, and it can be applied to any system.Nevertheless, proper calibration of the bio-optical model remains an essential factor for any site-specific application.The methodology presented here utilized relationships from a wide range of absorption and backscattering spectra that were presented by the IOCCG [3].This range was based on mathematical relationships that used phytoplankton concentration in ocean waters as a reference to estimate all absorption and backscattering components in the water column.Employing numerical modeling with bio-optical modeling validates the use of the presented approach for predictions of future scenarios related to changes in environmental and anthropogenic conditions.For example, environmental impacts due to global warming may cause alterations to seasonal cycles in temperature, precipitation and wind patterns.Such alterations can impact ocean circulation and stratification, which directly impacts the transport, distribution and composition of phytoplankton groups and consequently their optical properties.Similarly, anthropogenic effects due to increased aerosols can alter the incident light, which affects the photic depth and the wellbeing of the phytoplankton communities.

Discussion
A mathematical approach was presented to calculate irradiance in ocean waters (Case 1 waters).The methodology is not site-specific, and it can be applied to any system.Nevertheless, proper calibration of the bio-optical model remains an essential factor for any site-specific application.The methodology presented here utilized relationships from a wide range of absorption and backscattering spectra that were presented by the IOCCG [3].This range was based on mathematical relationships that used phytoplankton concentration in ocean waters as a reference to estimate all absorption and backscattering components in the water column.Employing numerical modeling with bio-optical modeling validates the use of the presented approach for predictions of future scenarios related to changes in environmental and anthropogenic conditions.For example, environmental impacts due to global warming may cause alterations to seasonal cycles in temperature, precipitation and wind patterns.Such alterations can impact ocean circulation and stratification, which directly impacts the transport, distribution and composition of phytoplankton groups and consequently their optical properties.Similarly, anthropogenic effects due to increased aerosols can alter the incident light, which affects the photic depth and the wellbeing of the phytoplankton communities.
It is argued here that the photic depth (at 1% irradiance) is an aggregate parameter for phytoplankton studies and that a more detailed representation of the "light quality" is more appropriate.As presented above, while the photic depth was 50 m, some wavelengths were extinct in much shallower depths (Figure 10).The resolved PAR distribution with depth (i.e., Figure 11) identifies light quality throughout the water column, which can provide useful information about the wellbeing of the phytoplankton communities at various depths.Such information can help in studies of primary productivity in the oceans.Most current numerical models calculate production (photosynthesis) based on total PAR only.In order for light quality to affect numerical models, they must compute the spectral photosynthesis of the competing dominant phytoplankton species and their adaptation to irradiance throughout the photic depth.An example of the complexity in this area is tackled for two coexisting phytoplankton species in [26].Partial numerical implementation was presented using λ = 490 nm as a proxy to the PAR [27].This topic needs more future studies to provide further guidance for implementing light quality and spectral photosynthesis in such models.
As indicated in the presented methodology, calibration is essential to reduce the level of uncertainty in the final irradiance predictions.Calibration of the bio-optical models eliminated the random processes that were incorporated in the original formulation [3].Constant values for the relevant mathematical parameters and coefficients were calibrated using field observations [24].The two physical parameters in Equation ( 9) can be site specific [3] and may have to be included in the calibration.Similarly, the spectral slopes (Equations ( 10) and ( 13)) have to be defined properly.Allowing such parameters to vary requires proper coupling between light and phytoplankton to account for the continually-changing parameterization (see the proposed numerical steps below).Validation of the presented methodology requires consistent high quality sets of data from various locations during various seasons.Figure 7 indicates that the observed downwelling irradiance failed to capture the attenuation within the top 12 m, which had the highest observed chl-a concentrations in the profile (Figure 4).Such inconsistency infers an unjustified increase in the observed downwelling irradiance.While physical effects (e.g., from density stratification) are beyond the scope of this work, they may have impacted the model and data shown.More data are definitely needed to improve the calibration and validation of the presented approach.
There are two major implications of the presented methodology: modeling implications and ecological implications.The methodology presented validates the coupling of bio-optical models with three-dimensional water quality models.The numerical implications of coupling phytoplankton to light in water quality models can execute the following general steps, which are recommended for relevant future work: 1.
At time t 2.
At location x,y 3.
For layer i 4.
The numerical model provides C x,y,z,t 5.
Calculate total irradiance, E i , by numerical integration over all λ (Equation ( 6)) 10.Move to the next layer down (i + 1) and repeat Steps 4-9 11.Move to next location and repeat Steps 3-10 12.For each layer, use E i in the calculation of C x,y,z,t+∆t , which fully couples light and phytoplankton at the next time level (t + ∆t).
The ecological implications encompass the predictive ability of phytoplankton biomass and its primary productivity in the oceans.Coupling light quality to water quality models facilitates understanding the relationship between light quality and phytoplankton.Light quality impacts the composition and health of phytoplankton communities in the water column, which, consequently, alter water clarity and irradiance levels.Water clarity is a major indicator used for many management decisions related to the health of estuarine, coastal and ocean waters.To the author's knowledge, such coupling is not present in most water quality models.However, recently, it was implemented in studies of primary productivity in the Pacific Ocean [27].
The presented approach is expected to work well in environments where phytoplankton particles are the major contributor to light attenuation through absorption by their chl-a pigment, their dead NAPs cells and their exudation of CDOM; as well as backscattering by the particulate phytoplankton cells and their dead NAPs.In other environments (e.g., with density stratification), the approach has to be augmented with additional considerations to account for such effects.Processes related to phytoplankton composition and dynamics (e.g., transport and vertical migration) should be covered by the water quality model.As the methodology is based on chl-a concentration, situations when phytoplankton cells and their chl-a content change may pose a challenge to this approach.

Conclusions
Predictive numerical models can use the presented methodology to couple light with phytoplnkton physiology throughout the photic depth.Not only light quantity but also its quality are essentisl for proper coupling.Nonetheless, numerical models have to accommodate the spectral photosynthesis of the competing dominant phytoplankton species and their adaptation to the spectral irradiance throughout the photic depth.The complexity of this coupling is still unresolved and needs more attention from both the numerical modeling community as well as the optical and bio-optical scientists.

Figure 1 .
Figure 1.Definition sketch of the relation between the numerical model vertical structure and its utilization to study light intensity throughout the water column.Concentration of chl-a (Cx,y,z,t) should be provided by the numerical model at location (x,y) at every layer depth (z) and at every time step (t) during the year.Ei−1 and Ei are the incident and departing downwelling irradiances through layer i (with thickness ℓi and extinction coefficient Ki) (Table1).

Figure 1 .
Figure 1.Definition sketch of the relation between the numerical model vertical structure and its utilization to study light intensity throughout the water column.Concentration of chl-a (C x,y,z,t) should be provided by the numerical model at location (x,y) at every layer depth (z) and at every time step (t) during the year.E i−1 and E i are the incident and departing downwelling irradiances through layer i (with thickness i and extinction coefficient K i ) (Table1).

Figure 2 .
Figure 2. Spectral distributions of absorption by pure water and the normalized absorption coefficients for phytoplankton pigment, NAPs and CDOM (full lines read on left axis); in addition to backscattering from phytoplankton, NAPs and water, as well as the derived spectral distribution function for the incident light (broken lines read on right axis).Refer to Table2for the definitions of normalized absorption and backscattering coefficients.

Figure 2 .
Figure 2. Spectral distributions of absorption by pure water and the normalized absorption coefficients for phytoplankton pigment, NAPs and CDOM (full lines read on left axis); in addition to backscattering from phytoplankton, NAPs and water, as well as the derived spectral distribution function for the incident light (broken lines read on right axis).Refer to Table2for the definitions of normalized absorption and backscattering coefficients.

2. 1 . 5 .
Light Backscattering by Phytoplankton, b p The calculations of backscattering for phytoplankton particles include the wavelength-dependent parameters for backscattering for the whole 400-700-nm spectrum, which are based on normalized values referenced to the wavelength λ = 550 nm [3].b p (λ) = b p b p (550) 1 are presented in [3].To eliminate the randomness, R 3 and P 3 are calibrated as presented in the Calibration and Validation Section.The original formulation in IOCCG [3] included the extra equation (b * p (λ) = b p (λ) − a p (λ)), which introduced superfluous error in b p when C x,y,z,t was zero.This equation is not included here.Examples of the values of the spectral distribution b + p (λ) = b p 550 λ n 1 for an arbitrary chl-a concentration are presented in Table

Figure 4 .
Figure 4. Summer data on 11 June 2013: (a) vertical distribution of ship-based measurements of chl-a concentrations and their vertical trend; (b) vertical distribution of glider-based average irradiance from 10:00 GMT-14:00 GMT; (c) spectral distribution of the near-surface average irradiance measurement.

Figure 4 .
Figure 4. Summer data on 11 June 2013: (a) vertical distribution of ship-based measurements of chl-a concentrations and their vertical trend; (b) vertical distribution of glider-based average irradiance from 10:00 GMT-14:00 GMT; (c) spectral distribution of the near-surface average irradiance measurement.

J 17 Figure 5 .
Figure 5. Fall data on 3 September 2013: (a) vertical distribution of ship-based measurements of chl-a concentrations and their vertical trend; (b) vertical distribution of glider-based average irradiance from 10:00 GMT-14:00 GMT; (c) spectral distribution of the near-surface average irradiance measurement.

Figure 5 .
Figure 5. Fall data on 3 September 2013: (a) vertical distribution of ship-based measurements of chl-a concentrations and their vertical trend; (b) vertical distribution of glider-based average irradiance from 10:00 GMT-14:00 GMT; (c) spectral distribution of the near-surface average irradiance measurement.

Figure 6 .
Figure 6.Calibration using fall data: (a) comparison between observed and predicted vertical irradiance profiles together with the expected irradiance in pure seawater; (b) correlation between observed and predicted irradiances.

Figure 7 .
Figure 7. Validation using summer data: (a) comparison between observed and predicted vertical irradiance profiles together with the expected irradiance in pure seawater; (b) correlation between observed and predicted irradiances.

Figure 6 .
Figure 6.Calibration using fall data: (a) comparison between observed and predicted vertical irradiance profiles together with the expected irradiance in pure seawater; (b) correlation between observed and predicted irradiances.

Figure 6 .
Figure 6.Calibration using fall data: (a) comparison between observed and predicted vertical irradiance profiles together with the expected irradiance in pure seawater; (b) correlation between observed and predicted irradiances.

Figure 7 .
Figure 7. Validation using summer data: (a) comparison between observed and predicted vertical irradiance profiles together with the expected irradiance in pure seawater; (b) correlation between observed and predicted irradiances.

Figure 7 .
Figure 7. Validation using summer data: (a) comparison between observed and predicted vertical irradiance profiles together with the expected irradiance in pure seawater; (b) correlation between observed and predicted irradiances.

Figure 8 .
Figure 8. Change of the total downwelling irradiance with depth on 3 September 2013.

Figure 9 .
Figure 9. Profiles and penetrations of various wavelengths, which illustrates the change in light quality with depth.

Figure 8 .
Figure 8. Change of the total downwelling irradiance with depth on 3 September 2013.

Figure 8 .
Figure 8. Change of the total downwelling irradiance with depth on 3 September 2013.

Figure 9 .
Figure 9. Profiles and penetrations of various wavelengths, which illustrates the change in light quality with depth.Figure 9. Profiles and penetrations of various wavelengths, which illustrates the change in light quality with depth.

Figure 9 .
Figure 9. Profiles and penetrations of various wavelengths, which illustrates the change in light quality with depth.Figure 9. Profiles and penetrations of various wavelengths, which illustrates the change in light quality with depth.

Figure 10 .
Figure 10.Penetration depth of the PAR on 3 September 2013.

Figure 11 .
Figure 11.Change of the spectral distribution of the PAR throughout the water column on 3 September 2013.

Figure 11 .
Figure 11.Change of the spectral distribution of the PAR throughout the water column on 3 September 2013.

Figure 11 .
Figure 11.Change of the spectral distribution of the PAR throughout the water column on 3 September 2013.

Table 1 .
Definitions of all abbreviations and symbols.

Table 2 .
Spectral distributions of absorption, a, backscattering, b, and the shape function, f, for the PAR (see the footnotes for details).
Column 1: wavelength counter, i; 9,10,16]tation location, (x-eastward, y-northward), and at the layer's vertical location z-below the water surface (m) and at the time, t, during the year.The coefficients 0.05 and 0.626 are based on observations from various regions including the North Atlantic, North Pacific, Gulf of Mexico, Mediterranean Sea, Arabian Sea and more.These coefficients can be site specific and can depend on λ [3,9,10,16].For convenience, the stated values of the two coefficients are used in the methodology presented here.The normalized phytoplankton absorption values, a c [13]), from Prieur and Sathyendranath[13](Table The range of the random n 2 values is −0.2-2.2, and its distribution is presented in [3].To eliminate the randomness, R 4 and P 4 are calibrated as presented in the Calibration and Validation Section.Examples of the values of the spectral distribution b + s (λ) = b s 550 λ n 2 for an arbitrary chl-a concentration are presented in Table