Capturing optically important constituents and properties in a marine biogeochemical and ecosystem model

Abstract. We present a numerical model of the ocean that couples a three-stream radiative transfer component with a marine biogeochemical–ecosystem component in a dynamic three-dimensional physical framework. The radiative transfer component resolves the penetration of spectral irradiance as it is absorbed and scattered within the water column. We explicitly include the effect of several optically important water constituents (different phytoplankton functional types; detrital particles; and coloured dissolved organic matter, CDOM). The model is evaluated against in situ-observed and satellite-derived products. In particular we compare to concurrently measured biogeochemical, ecosystem, and optical data along a meridional transect of the Atlantic Ocean. The simulation captures the patterns and magnitudes of these data, and estimates surface upwelling irradiance analogous to that observed by ocean colour satellite instruments. We find that incorporating the different optically important constituents explicitly and including spectral irradiance was crucial to capture the variability in the depth of the subsurface chlorophyll a (Chl a) maximum. We conduct a series of sensitivity experiments to demonstrate, globally, the relative importance of each of the water constituents, as well as the crucial feedbacks between the light field, the relative fitness of phytoplankton types, and the biogeochemistry of the ocean. CDOM has proportionally more importance at attenuating light at short wavelengths and in more productive waters, phytoplankton absorption is relatively more important at the subsurface Chl a maximum, and water molecules have the greatest contribution when concentrations of other constituents are low, such as in the oligotrophic gyres. Scattering had less effect on attenuation, but since it is important for the amount and type of upwelling irradiance, it is crucial for setting sea surface reflectance. Strikingly, sensitivity experiments in which absorption by any of the optical constituents was increased led to a decrease in the size of the oligotrophic regions of the subtropical gyres: lateral nutrient supplies were enhanced as a result of decreasing high-latitude productivity. This new model that captures bio-optical feedbacks will be important for improving our understanding of the role of light and optical constituents on ocean biogeochemistry, especially in a changing environment. Further, resolving surface upwelling irradiance will make it easier to connect to satellite-derived products in the future.

. Spectra for (a) absorption and scattering by water molecules (a w , b w , m −1 ); (b) particle-specific absorption and scattering by detritus (a part det , b part det , m 2 particle −1 ); (c) CDOM-specific absorption by CDOM (a CDOM cdom , m 2 (mmol P) −1 ); (d) Chl a-specific total absorption by phytoplankton (a chl phyj , m 2 (mg Chl a) −1 ); (e) Chl a-specific absorption by photosynthetic pigments (a chl psj , m 2 (mg Chl a) −1 ); and (f) biomassspecific scattering by phytoplankton (b C phyj , m 2 (mg C) −1 ). Details on data sources are included in the main text and Appendix Sect. C. The black line in (d-f) is the mean of the coloured lines (i.e. the mean spectrum). Spectra are shown here with 1 nm resolution for clarity; the model uses the average over the 25 nm bands (vertical grey lines).
ter. Attenuation of light within the water column is an interaction of absorption and scattering by "optically important constituents", including water molecules, detrital matter, coloured dissolved organic matter (CDOM), and the phytoplankton themselves.
Phytoplankton absorb light in the visible spectrum (400 and 700 nm). The optical constituents attenuate these wavelengths differently. For instance, water molecules absorb very strongly in the longer wavelengths ( Fig. 1a), while detrital matter and CDOM absorb more in the shorter wavelengths (Fig. 1b,c). Thus the spectrum of light at any location is a complex function of the combination of different optical constituents in the overlying water. Previous studies have highlighted the importance of resolving the spectral light field (e.g. Fujii et al., 2007;Kettle and Merchant, 2009), especially as different species of phytoplankton have different light absorption spectra (e.g. Stramksi et al., 2001;Sathyendranath and Platt, 2007). This difference in efficiency of light absorption by phytoplankton is important for their relative fitness and biogeography (Bidigare et al., 1990a;Huisman and Weissing, 1995;Moore and Chisholm, 1999;Stomp et al., 2004;Hickman et al., 2010).
Much is known about the optics of water (e.g. Pope and Fry, 1997;Smith and Baker, 1981;Morel, 1974;Zhang and Hu, 2009;Kirk, 1994). Although much is known about the distributions of CDOM (Nelson and Siegel, 2013), detritus (Loisel, 2002), and phytoplankton (IOCCG report 15, 2014) it remains unclear how their distributions feed back to phytoplankton community structure and biogeochemistry. Numerical models provide useful tools to explore these interactions, but to do so requires an appropriately detailed description of the photosynthetically available radiation (PAR).
Several recent models resolve the light spectrum and some of the absorption and scattering properties of different constituents (e.g. Mobley et al., 2009;Fujii et al., 2007;Gregg and Casey, 2007;Bisset et al., 1999). Such models include fully coupled radiative transfer, but differ in the levels of simplification for computational efficiency (e.g. Fujii et al., 2007;Gregg and Casey, 2007) and differ in which and how they treat the different water constituents. For instance, CDOM is treated as uniform in Fujii et al. (2007), and linked to chlorophyll a (Chl a) in Gregg and Casey (2007). Fujii et al. (2007) suggested that including explicit optics in an ecosystem model allowed a more accurate subsurface light field as well as additional constraints on model parameters. Several additional studies have demonstrated the value of adding optics to biogeochemical models (e.g. Babin et al., 1993;Sathyendranath and Platt, 2007;Kettle and Merchant, 2008).
In Sect. 2 we introduce an updated version of the MIT biogeochemistry and ecosystem model (Follows et al., 2007;Dutkiewicz et al., 2012) with a radiative transfer component as well as the explicit treatment of several optical constituents (water molecules, detrital matter, CDOM, S. Dutkiewicz et al.: Modelling optical properties 4449 and a community of optically distinct phytoplankton types). Specifically, each constituent is treated independently. The fully coupled radiative transfer allows us to calculate spectral surface upwelling irradiance; a product similar to that measured by ocean colour satellites. We show results from this new coupled model where the light field is a dynamic function of the different optical constituents and evaluate against several data sets (Sect. 3). In particular, we use a comprehensive data set from an Atlantic Meridional Transect cruise which includes detailed concurrent optical, biogeochemical, and ecosystem observations between the UK and South Africa in September/October of 2004 (AMT-15). Some of the observations are published here for the first time. The data set is ideal for evaluating how our model captures the amount and nature of the light that penetrates the water column across basin scale along with the relevant ecological properties.
We perform a number of sensitivity experiments that explore the value of the additional model complexity (Sect. 4), the role of each of the water constituents (Sect. 5), and their relative importance. The model allows us to investigate how changes to any constituent feed back to the system, impacting phytoplankton biogeography, biogeochemistry, and surface irradiance reflectance.

Model description
The biogeochemical-ecosystem model resolves the cycling of carbon, phosphorus, nitrogen, silica, iron, and oxygen through inorganic, living, dissolved, and particulate organic phases as discussed in Follows et al. (2007), Dutkiewicz et al. (2009Dutkiewicz et al. ( , 2012, and Hickman et al. (2010). The biogeochemical and biological tracers are transported and mixed by a the MIT general circulation model (MITgcm) (Marshall et al., 1997). The physical framework is flexible, but here we employ a global configuration which is constrained to be consistent with altimetric and hydrographic observations (the ECCO-GODAE state estimates; Wunsch and Heimbach, 2007). This three-dimensional configuration has 1 • × 1 • horizontal resolution and 23 levels ranging from 10 m in the surface to 500 m at depth. These physical fields have been used in many previous biogeochemical-ecosystem studies (e.g. Follows et al., 2007;Dutkiewicz et al., 2009Dutkiewicz et al., , 2012Prowe et al., 2012).
Similar to several of these previous studies, we resolve several phytoplankton types, P j , as well as two simple grazers, Z k . The biogeochemical and biological tracers interact through the formation, transformation, and remineralization of organic matter. Excretion and mortality transfer living organic material into sinking particulate and dissolved organic detritus, which are respired back to inorganic form. Aeolian iron fluxes to the ocean surface are provided by Luo et al. (2008). We provide complete model equations, description and parameter values in Appendix Sect. A and Tables 1 to 6. Here we focus on the relevant new features: in particular an explicit radiative transfer component that allows us to consider absorption and scattering of light spectrally and with attention to each of the relevant optical constituents.

Radiative transfer model
Irradiance just below the surface of the ocean is provided by the Ocean-Atmosphere Spectral Irradiance Model (OASIM) (Gregg and Casey, 2009)  ). OASIM includes the impact of clouds, water vapour, and aerosols in the atmosphere and surface roughness and reflectance at the ocean-atmosphere interface. Irradiances are provided averaged in 25 nm wavebands from 400 to 700 nm. The two downward light streams (direct and diffuse, E d , E s ) in each waveband are followed through the water column. Irradiance is attenuated by absorption (a),and scattering (b), which includes both forward (b f ) and backward (b b ) components. Scattering diverts irradiance from the direct and diffuse beams and partitions it between the downward diffuse and upwelling stream (E u ).
We parameterize this "three-stream" irradiance model following Aas (1987), Ackleson et al. (1994), and Gregg (2002). The model is described by the simultaneous equations for the light streams in each waveband (λ) with depth (z): where r s , r u , and r d are the effective scattering coefficients that are normalized by backward scattering coefficients; υ d , υ s , and υ u are the average cosines (definition in Appendix Sect. B); and the radiance is separated into a direct beam and a diffuse component. This set of equations can be simplified following Aas (1987) by approximating r s , r u , r d , υ s , and υ u with constant values (see Appendix Sect. B). With these assumptions, the set of equations can be reduced to a tri-diagonal system. In contrast to Aas (1987), Ackleson et al. (1994), and Gregg (2002) we solve E d (λ), E s (λ), and E u (λ) directly at each location and at each depth using Gaussian elimination.
We calculate total scalar irradiance, E 0 (λ), in each waveband at each location and layer (averaged, multiplicatively, between the top and bottom) by scaling the irradiance by the This is the light available to the phytoplankton. We note that the radiative transfer component is a simplification from a full radiance model and, in particular, does not resolve the angular distribution of light or angular dependence of scattering. These assumptions have been shown to be small in terms of the needs for ecosystem models (Mobley et al., 2009). Though not a full radiative transfer model, our three-stream treatment does provide the relevant output for our objectives: the spectral light available for photosynthesis and an upwelling component that at the sea surface is similar to that seen by a satellite.

Surface reflectance
Since the model resolves an upwelling stream of irradiance, we can calculate a surface reflectance (unitless): where E below u (λ)| k=0 is upwelling irradiance just below the surface and E below d o (λ) + E below s o (λ) represents the downward (direct and diffuse) irradiance just below the surface as provided by OASIM.
To compare to remotely sensed reflectance (R RS ) we convert between model subsurface reflectance and the slant upward radiance seen by satellite by using a bidirectional function Q: The bidirectional function Q has values between 3.5 and 5 sr depending on many variables, including inherent optical properties of the water, wavelength, and solar zenith angles Voss et al., 2007). For simplicity here we assume that Q = 4 sr. Model R RS is therefore analogous but not exactly the same as that measured by satellite. R RS has units of sr −1 .

Treatment of water constituents
Attenuation of irradiance results from absorption by water molecules (a w ), phytoplankton (a phy ), detrital particles (a det ), and CDOM (a cdom ) and from scattering by water molecules (b w ), phytoplankton (b phy ), and detrital particles (b det ). The absorption (a), total scattering (b), and backward scattering (b b ) (all with units of m −1 ) are represented as In the model we use absorption and scattering coefficients ( Fig. 1) averaged over 25 nm bands to match the irradiance input from a variety of sources, detailed below.

Water molecules
We assume absorption by water molecules (a w , b w , b bw ) to follow the spectrum of Pope and Fry 1997. Scattering is taken from Smith and Baker (1981) and Morel (1974), and backscattering from Morel (1974) and Morel et al. (2007). The spectra for these are shown in Fig. 1a.

Detrital matter
The model uses the absorption and scattering spectra for detrital matter (Fig. 1b) from Stramski et al. (2001). These spectra were deduced by assuming an assemblage of particles with size distribution described by a power function with slope of −4, and the values are given in terms of absorption or scattering per particle (Stramski et al., 2001). Thus we introduce the coefficient p part to convert the model particulate organic carbon (POC) to number of particles, making the crude assumption that the size distribution of particles is uniform everywhere. The absorption and scattering by particles is described as Here we use the convention that the superscript on the a, b, and b b terms refers to the normalization variable, here particle concentration. Units of a part det (λ), b part det (λ), and b part bdet (λ) are m 2 particle −1 .
We note that in the optical community the term "nonalgal particles", or NAP, is frequently used for any nonphytoplankton particles. In this paper we specifically use the term "detritus" instead, as we link to the non-living organic matter pool and do not explicitly resolve other non-algal particles such as viruses and heterotrophic bacteria.

Coloured dissolved organic matter
CDOM absorbs highly in the short wavelengths and absorption decreases exponentially with increasing wavelength (Kitidis et al., 2006;Nelson and Siegel, 2013). CDOM is not usually explicitly resolved in marine ecosystem models (exceptions are Xiu andChai, 2014, andBissett et al., 1999). Here we have resolved an explicit CDOM-like tracer (denoted "CDOM") similar to Bissett et al. (1999). The model CDOM has units of concentration (mmol C m −3 ), and is assumed to have a source that is a fraction (f cdom ) of DOM production, to have a long remineralization timescale (d cdom ), and to be bleached under high light conditions. The bleaching is parameterized to reach a maximum rate, ι cdom , when PAR is above I cdom , and linearly decreases at lower PAR. The sources and sinks of this CDOM-like tracer are there-  Elemental ratios where S DOM S represents the sources of DOM (see Appendix Sect. A) and γ T is the temperature function affecting biological rates. We parameterize a cdom (λ) as function of "CDOM" such that and where a CDOM cdom (λ) is the concentration-specific absorption of the CDOM-like tracer (Fig. 1c). The value for the spectral slope, s cdom , is taken from the literature (Kitidis et al., 2006), and c cdom (λ o ) is the CDOM-specific absorption at reference waveband, λ o . Although CDOM is also strongly linked to terrestrial matter, we do not provide any land sources at present. We discuss the sensitivity of the function and parameters, and compare to previous model parameterizations in Sect. 5.

Phytoplankton
The absorption and scattering by phytoplankton is the net effect of each phytoplankton type resolved in our model, j : The Chl a-specific absorption spectra a chl phy j (λ) have units of m 2 (mg Chl a) −1 . The scattering (b C phy j (λ)) and backscattering (b C bphy j (λ)) coefficients are assumed to be functions of phytoplankton biomass (e.g. Martinez-Vincent et al., 2013) and has units m 2 (mol C) −1 . These spectra are specific to each of the phytoplankton types j (Fig. 1d-f) as taken from the literature (see discussion in Sect. 2.5 and Appendix Sect. C). M Cj is the C : P ratio in each phytoplankton type (the base "currency" of the plankton in the equations is phosphorus; see further details in Appendix Sect. A).

Phytoplankton growth
Phytoplankton growth is modelled as a function of temperature, irradiance, and nutrients as in Hickman et al. (2010) following Geider et al. (1998). The growth rate is equal to the carbon-specific photosynthesis rate: where P C mj is the light-saturated photosynthesis rate that is a function of temperature and nutrient limitation (see Appendix Sect. A), and θ j is the ratio of Chl a to C within each phytoplankton j (discussed further below). Furthermore, is the scalar irradiance absorbed by each phytoplankton, j , multiplied by φ max j , the maximum quantum yield of carbon fixation. Ej is thus equivalent to the spectrally resolved product of irradiance and the initial slope of the Chl a normalized photosynthesis versus irradiance curve and has units of mmol C (mg Chl a) −1 d −1 . a chl psj (λ) is the Chl aspecific photosynthetic absorption spectra in each waveband λ (Fig. 1e), and E 0 (λ) comes from the radiative transfer code (see Eq. 4).
Since some pigments are photoprotective, phytoplankton do not use all the light that they absorb for photosynthesis. Similar to Hickman et al. (2010) and Bisset et al. (1999), the total absorption spectra are therefore greater than the photosynthetic absorption spectra, a chl phyj > a chl psj (Fig. 1d, e). (See Max. grazing rate g max j k d −1 DOM-POM partitioning ϕ g ij k unitless ϕ mz ik unitless Mortality at 30 • C m zk d −1 m z2k d −1 m 3 (mmol P) −1 Grazing efficiency ζ j k unitless Grazing half-saturation κ pk mmol P m −3 discussion in Sect. 2.5.) We also allow for photoinhibition, as in Hickman et al. (2010), such that P C mj reduces above a critical value at high light (see Appendix Sect. A).
Cell size governs many traits. Smaller phytoplankton have lower nutrient half-saturation constants and sink more slowly. The maximum growth rates are guided by observations, diatoms having the highest rates and Prochlorococcus having the lowest (see e.g. Irwin et al., 2006). The parameter values are within ranges found in the literature and previous ecosystem model Ward et al., 2013;Monteiro et al., 2010).
In this model we treat the phytoplankton light absorption and scattering explicitly (Sect. 2.3.4). The Chl a-specific absorption spectra a chl phyj (λ) (units, m 2 (mg Chl a) −1 ) vary between functional group (Fig. 1d). These spectra were obtained from representative phytoplankton types in cultures grown at similar growth irradiance (see references in Appendix Sect. C). The spectra capture differences in pigment composition and other taxon-specific differences, including the "package effect" (Berner et al., 1989). For instance, the larger diatom has a flatter spectrum than the smaller phytoplankton (e.g. Prochlorococcus). Total light-scattering spectra (b C phyj , Fig. 1f) were also obtained from representative phytoplankton types in culture, as were the ratios of backscatter to total scatter for each phytoplankton (b C bphyj , units m 2 (mol C) −1 ) (Stramski et al., 2001;Subramaniam et al., 1999). Spectra for absorption by photosynthetic pigments (a chl psj , Fig. 1e) were derived using the pigment reconstruction technique (following Hickman et al., 2010;Babin et al., 1996). Light absorption spectra were reconstructed by scaling  the weight-specific absorption coefficients for Chl a, b, and c; photosynthetic carotenoids and non-photosynthetic carotenoids; phycoerythrobilin; and phycourobilin-rich phycoerythrins (Bidigare et al., 1990b) to obtain the lowest sum of residuals between reconstructed and observed spectra. a chl psj was then calculated by adjusting the measured a chl phyj by the spectral ratio of the reconstructed spectra with and without non-photosynthetic pigments (Hickman et al., 2010).
We parameterize all phytoplankton to have the same maximum quantum yield of carbon fixation (φ max j , units mol C fixed per moles photons) and all but diatoms to have the same maximum Chl a : C (θ max j , units mg Chl a (mmol C) −1 ) (MacIntyre et al., 2002). We parameterize low-light Prochlorococcus as being photoinhibited, as this is a distinct feature of the difference between highand low-light strains (Moore and Chisholm, 1999;Hickman et al., 2010).
We resolve two zooplankton classes (large and small) that graze on the phytoplankton using a Holling III scheme (Holling, 1959). The large class preys preferentially on the diatoms, coccolithophores, and Trichodesmium, while the smaller class preys preferentially on the smaller phytoplankton. We additionally parameterize diatoms and coccolithophores (hard shells) and Trichodesmium (toxicity) as having lower palatability. Zooplankton grazing parameters are similar to those used in Prowe et al. (2012), which were determined from a mechanistic model of zooplankton feeding (see Table 6).

Enhancements and limitations of optics component
The inclusion of radiative transfer and spectral light, as well as capturing several important optical constituents, is a significant development in the model. However, this version of the model is not without limitations. One major, though currently necessary, simplification is to assume constant absorption and scattering spectra ( Fig. 1) for each constituent. For instance, absorption spectra for phytoplankton types do in reality change based on shifts in Chl a : C (e.g. MacIntyre et al., 2002;Morel et al., 1993Morel et al., , 1995 as well as changes in ratios of photoprotective to photosynthesis pigments as a result of light, temperature, and nutrient stress (e.g. Stramski et al., 2002). However, these changes are likely to be small compared to the differences already captured by the representative spectra and photoacclimation component, and there are not, as yet, enough systematic observations of all of these alterations to constrain model parameterizations. Additionally, the CDOM absorption spectrum has been observed to alter regionally (e.g. Kitidis et al., 2006;Twardowski et al., 2004;Bricaud et al., 2010), though as yet we feel it is premature to attempt to capture this variability in the model parameterizations.
Scattering, particularly by detrital particles, remains the least well developed aspect of the model. In particular, we neglect variations in detrital particle size distributions, which are likely to be important (Stramski et al., 2001). Additionally, the spectrum for b part det that we use (Stramksi et al., 2001, Fig. 1b) makes the assumption of homogeneous spheres. However it is likely that differences in shapes and internal structure of the particles will be important for altering the spectral shape (Stramski et al., 2004). We also do not take into account inelastic scattering, which may be important for blue and green light in oligotrophic regions (e.g. Ge et al., 1993).
We additionally currently neglect other potentially important optical constituents such as minerals (e.g. Stramski et al., 2001), particulate inorganic carbon (e.g. Balch and Itgoff, 2009), colloids and bubbles (e.g. Stramski et al., 2004), and non-photosynthetic organisms including zooplankton, bacteria (e.g. Morel and Ahn, 1991), and viruses (e.g. Stramski et al., 2001). We felt that these are, as yet, not well enough constrained to include explicitly in the model. The limitations list above should not, however, detract from the major enhancement to the model, and our assumptions are similar to those of other models (e.g. Fujii et al., 2007;Gregg and Casey, 2007). This new model provides a unique platform to examine global implication of optical properties to the phytoplankton ecosystem, feedbacks to the biogeochemistry, and links to satellite data that are not possible with limited observational data. Here we first validate the model in a standard "default" configuration. We then provide a series of studies exploring the significance of each of the optical constituents and our parameterization. Several studies in progress build on these results.

Default simulation and validation
We initialize the macronutrient fields (nitrate, phosphate, and silicic acid) from World Ocean Atlas (Garcia et al., 2006) climatologies and the iron from previous model output. We also use previous model output to provide distribution of the ammonium, nitrite, and dissolved and particulate matter. The total phytoplankton biomass is initialized from previous model output, divided equally between groups, except for the diazotrophs, which are initialized at a much lower value so as not to flood the system with new nitrogen in the first few time steps. Zooplankton are similarly initialized with equal distribution in both groups.
The model time step is 3 h. We tested this against smaller time steps with almost identical results. We run the simulation forward for 10 years with a repeating generic "year" from the physical ECCO-GODAE products (Wunsch and Heimbach, 2007). Model results shown in this section are from the last year of the simulation. The phytoplankton establish a repeating pattern after about 3 years such that we can assume a "quasi-steady state" by year 10. A slow drift as deep water nutrient distributions adjust does not significantly change the results over the remaining time period. We evaluate the model results against a range of in situ observations and satellite-derived products. In particular, we focus on the unique data set including biogeochemical, ecological and (some previously unpublished) optical properties that were obtained as part of the AMT-15 cruise. Though there are other AMT cruises that include some similar and/or different combinations of optical data (e.g. AMT-19, Dall'Olmo et al., 2012, Martinez-Vicente et al., 2013, we chose to look at only a single transect for clarity. In particular, the combination of data on spectral irradiance penetration, a CDOM , and light absorption by phytoplankton was of particular use in model validation.

Atlantic Meridional Transect
The model broadly reproduces the horizontal gradients at the surface but, importantly, also captures the subsurface Chl a maximum ( Fig. 3a, b), and in particular its deepening in the subtropical gyres, especially in the South Atlantic. The model captures the depth of the nitricline across the transect (Fig. 3c, d), especially the deep section (200 m) in the South Atlantic gyres. The model does not adequately resolve the North Atlantic upwelling (a resolution issue in the physical model), and nitrate and Chl a are too low in this region. Additionally, the physical model has too strong upwelling just south of the Equator, leading to nitrate and Chl a being too high.
The model also captures observed variability in a cdom along the AMT-15 transect: low in the surface waters where CDOM is quickly bleached, and higher in deeper waters where CDOM accumulates. Values and regional patterns compare well between model and observations (Fig. 3e, f), except just south of the Equator, where Chl a and nutrient supply are also too high (as discussed above). Absorption by phytoplankton ( Fig. 3g) was only measured at the surface and the subsurface Chl a maximum. The model captures the higher value near the subsurface Chl a maximum (Fig. 3h).
We have used the AMT-15 measured downwelling irradiance and upwelling zenith radiance together with the inversemodelling procedure of Boynton (1997, 1998) to estimate the total absorption and total backscattering in several wavelengths (Fig. 4a, c, e, g). We discuss this inversion further in Appendix Sect. D. There is a large degree of uncertainty in this inversion process and additional noisiness provides several spurious high/low values that are not realistic. Given this caveat, we find that the model qualitatively captures (Fig. 4b, d, f, h) the magnitudes and the pattern of higher absorption/lower scattering at the higher wavebands. Since the model realistically captures much of the variability in optical constituents it also accurately resolves the penetration of light through the water column (Fig. 5) as found in the AMT-15 data. We compare the depth of the 1 % light level: the depth where the downwelling irradiance in each waveband is 1 % of the surface value (E below . We find that the shortest wavebands (e.g. purple line and symbols in Fig. 5) reach deepest in the South Atlantic gyre, where concentrations of the optical constituents are lowest and less deep than medium wavebands (e.g. light-and darkblue lines) in more equatorial regions. The penetration of blue wavebands leads to the very deep subsurface Chl a maximum and drawdown of nutrients at depth as observed in the AMT-15 transect and in the model. The 1 % depths are too deep in the North Atlantic upwelling region, since we do not capture this feature in the physics.
The model captures intricate patterns of absorption and scattering that develop from the interplay of different optical constituents and suggests the importance of treating each constituent separately for reproducing the in situ light field. We explore this further in Sect. 5.

Global results
That the model captures much of the Chl a, nutrient, and optical properties at basin scale and with depth as observed during AMT-15 is very encouraging. The model also captures many of the global features in Chl a (derived from MODIS satellite), primary production (derived using Behrenfeld and Falkowski, 1997), and macronutrients (from the World Ocean Atlas; Garcia et al., 2006), though with notable biases (Fig. 6). The broad-scale features of high nutrients, high Chl a, and high productivity in the high latitudes and equatorial regions as well as low nutrients and low Chl a in the subtropical gyres are resolved. We do not, however, capture coastal features as the physical model is too coarse to resolve the important mesoscale processes. This is also true in frontal zones (such as the western boundary currents) where primary production is too low.
Relative to the composite of iron data (Tagliabue et al., 2012), we also capture high iron in the Atlantic Ocean and lower iron over much of the Pacific (Fig. 6j, k). However, iron may be too low in the tropical South Pacific and Pacific equatorial regions. Here the model aeolian dust supply (based on Luo et al., 2008) may be too low; however the physical model also does not adequately resolve equatorial undercurrents, which are likely responsible for supplying sedimentary iron to this region (Radic et al., 2011;Slemons et al., 2009). Since iron limitation is too strong in this region, productivity is too low and nitrate too high. The model also overestimates Chl a in the Southern Ocean and other high latitudes relative to the satellite product. However, the satellite Chl a algorithm has a factor of 2 range error (Campbell et al., 2002) and is especially problematic in the Southern Ocean (Szeto et al., 2011).
We find that the spatial standard deviation (between 0.85 and 1.15) and correlation (greater than 0.9) of the model vs. observed nutrients are encouraging (Fig. 7). Though we capture much of the spatial variability in the Chl a, the corre-lations with satellite-derived products are not as good. The primary production is universally too low and too uniform relative to the satellite-derived product. However, we note that the satellite products of Chl a and primary production have large error margins associated with them that are not spatially homogeneous (Szeto et al., 2011).
The model ecosystem has distinctive seasonal cycles (Fig. 8) that mostly match the observed satellite-derived and in situ Chl a at nine time series sites (locations shown in Fig. 2) collected as part of JGOFS (Kleypas and Doney, 2001). In many locations the model overestimates the satellite-derived peak of the bloom (consistent with annual mean Chl a being too high) but captures the non-bloom values more accurately. However, the in situ data broadly encompass the model values. We also capture the satellitederived timing of the spring bloom, though notably do not get correct blooms at Station P, Kerfix, NABE.
A unique feature of this model is irradiance reflectance output, which we have converted to remotely sensed reflectance (R RS ) using a fixed bidirectional function Q (see Sect. 2.2). We compare this model output to MODIS remotely sensed reflectance, R RS (λ). Despite the mismatch in wavelength and bandwidth and the oversimplification of a fixed Q, the model qualitatively captures the pattern of high reflectance in the subtropics relative to the higher productivity regions in low wavebands and the opposite pattern in higher wavebands. These initial results suggests that the model framework will be a useful laboratory for exploring satellite-like semi-analytical inversion algorithms (e.g. IOCCG report 5, 2006).

Phytoplankton biogeography
Eight of the nine phytoplankton functional groups that we resolve have distinct biogeography (Fig. 10). This biogeography encompasses both horizontal and vertical patterns of phytoplankton biomass. The large eukaryote group does not survive in this model as it was given no specific trade-off. It was large (low nutrient affinity) and had a low growth rate (typical of dinoflagellates).
We compare simulated biomass of the picophytoplankton to observations from AMT-15 (Fig. 11). AMT-15 cell counts were measured by analytical flow cytometry following methods of Heywood et al. (2006) and converted to biomass using constant factors (Zubkov et al., 1998) for comparison purposes. The model captures the smallest autotrophs, Prochlorococcus, as having significant abundances through the subtropics and tropic; Synechococcus was more abundant at the northern poleward fringe of the subtropics, and picoeukaryotes were more ubiquitous and more dominant in the subsurface Chl a maximum. In the 20 to 5 • S region the model nutrient source is too high (discussed above) and Synechococcus analogues unrealistically dominate instead in the model. The model distribution of large phytoplankton biomass (e.g. diatoms, coccolithophores) compared well to observations made along other AMT cruises Cermeño et al., 2008).
The MAREDAT (MARine Ecosystem DATa; Buitenhuis et al., 2013) compilation provides a comprehensive, though still sparse, climatological distribution of several plankton functional groups. Here we re-grid the MAREDAT compilation onto a 5 • grid with all observations between 0 and 50 m averaged together and compare this to the model output (Fig. 12). For the model results we sum the Prochlorococcus, Synechococcus, and picoeukaryote groups to compare to the observations of picophytoplankton. The model captures the ubiquitous nature of the picophytoplankton, the lack of coccolithophores in the subtropical gyres and polar extent of the Southern Ocean, the high diatom values in the high latitudes and in the equatorial upwelling regions, and the low abundance (or lack) of diazotrophs in the southern Pacific gyres. However, model coccolithophore biomass is in general too high and diazotroph biomass has a peak too far south in the North Atlantic. Though the MAREDAT compilation includes micro-, meso-, and macrozooplankton, data on the former and the latter are very sparse. Since we do not have direct analogues in the model, we show here only the mesozooplankton biomass observations (Fig. 12i). The model captures the patterns of high and low values of zooplankton biomass.
Given the sparsity of in situ measurements of phytoplankton types, it is natural to attempt to capture aspects of bio- geography from space (IOCCG report 15, 2014; IOCCG report 9, 2009). Here we compare the model output to the PHYSAT product (Alvain et al., 2008), which empirically relates optical properties to specific (probably dominant) phytoplankton types (Fig. 13a, c) for January and July and compare to model dominant types (Fig. 13). In both the model and PHYSAT we find that cyanobacteria dominate the tropics and subtropics, diatoms play a substantial role in the summer biomass, and a combination of coccolithophores and picoeukaryotes dominate in the mid-latitudes.
These global "observations" contain many uncertainties stemming mainly from the scarcity of in situ data, but the model does not disagree with their findings. The model captures key patterns of observed optical and ecological properties. It provides a tool to explore aspects of the ocean biogeochemistry and ecosystem that are not possible with models that do not explicitly resolve radiative transfer, spectral irradiance, and resolution of the different water optical properties. In the next section we explore the role of the various water constituents on the irradiance spectrum and how they impact biogeochemistry and ecosystem structures.

Sensitivity experiments: value of added model complexity.
We conduct two sensitivity experiments to highlight the importance of the extra levels of complexity of this new version of the model. In the first experiment (designated EXP-V0) the biogeochemistry and ecosystem are the same as in the default experiment described above (designated EXP0) but there is only a single band of irradiation (400-700 nm, summed over the original 25 nm, so that total PAR is conserved); attenuation (c tot ) of PAR is a function of only absorption by water molecules and Chl a summed over all phytoplankton types: c tot = a wo + a chlo Chl a tot , where a wo = 0.04 m −1 and a chlo = 0.04 m 2 (mg Chl a) −1 . There is no explicit account taken for optical role of CDOM or detritus (though the value chosen for a chlo does implicitly include their role). Similar parameterizations have been used in previous versions of our model (e.g. Dutkiewicz et al., 2014) and are also common in many other biogeochemical models.
The results from EXP-V0 (Fig. 14a) reveal a much more latitudinally uniform penetration of light. In particular the subsurface Chl a maximum in the subtropical gyre is too shallow relative to the default experiment (EXP0, Fig. 14c) and observations (Fig. 3a).
In experiment EXP-V1 we include all the optical constituents explicitly (as in EXP0), though with only a single band of PAR (as in EXP-V0). We assume the absorption and scattering coefficients for 500 nm in this experiment. This experiment (Fig. 14b) reveals a substantially more realistic variability in the depth of the subsurface Chl a maximum and penetration of PAR. The addition of spectral light leads to even deeper penetration of PAR in the subtropical gyres (compare to default experiment, EXP0, Fig. 14c): deepest penetrating light is in the blue/green range, such that an average absorption across one waveband will not capture these differences.
These sensitivity experiments suggest that explicitly capturing regional changes in all optical constituents is essential for the realistic variations in the depth of light penetration. Resolving the light spectrum further enhances the realism of the results. The addition of the radiative transfer code is essential for obtaining upwelling irradiance that can link to satellite products.

Sensitivity experiments: role of optical constituents
Optical constituents play varying roles in their effect on irradiance attenuation (absorption and scattering). These roles have long been a topic of interest; however many studies have included only limited observations and been of highly localized in character (e.g. Jerlov, 1953;Chang and Dickey, 1999) but have nonetheless recognized that they vary regionally (e.g. Barnard et al., 1998;Simeon et al., 2003). Targeted cruises have also provided larger-scale observations indicating a wide range of values for each constituent and altering importance in different regions (e.g. BIOSOPE; Bricaud et al., 2010). Additionally, several attempts have been made to construct algorithms to determine the relative contributions from more easily measured quantities, including those from satellites (e.g. Maritorena et al., 2002;Lee et al., 2002Lee et al., , 2007Ciotti and Bricaud, 2006;Werdell et al., 2013;Zheng and Stramski, 2013). Our model provides a unique global three-dimensional perspective. Here our results focus on an (extended) AMT transect (Figs. 15 and 16); however, they are also consistent with observations in other regions (e.g. Bricaud et al., 2010).
Absorption by water molecules is most important at longer wavebands (Pope and Fry, 1997) but still has an impact at shorter wavebands (Fig. 15a, b, i, j). It is relatively more important in lower productive waters (e.g. South Atlantic gyre) because the concentrations of other constituents are relatively low. Absorption by detrital matter plays a role, especially near the 1 % depth in highly productive regions and at shorter wavebands (Fig. 15c, d, i, j) as suggested by observations (e.g. Jerlov, 1953). Absorption by phytoplankton plays a significant role where Chl a is highest (e.g. the subsurface Chl a maximum, as found in observations; e.g. Chang and Dickey, 1999) at wavelengths less than 550 nm, and little role at longer wavelengths (Fig. 15g, h, i, j; see also Fig. 1). Absorption by CDOM at short wavebands is important (as seen in observations; e.g. Jerlov, 1953) in most regions, particularly where productivity is high where it is the dominant absorber. It also has, relative to other constituents, a large role at depth (as seen in observations; e.g. Simeon et al., 2003;Bricaud et al., 2010;Nelson and Siegel, 2013). At long wavebands, CDOM plays very little role. Scattering by   phytoplankton is most important at shallower depths, while scattering by detrital matter is dominant deeper at all wavelengths (Fig. 16).
We perform a series of sensitivity experiments to explore the role of each constituent in setting the irradiance field in the ocean and on surface reflectance, and to see how changes to these constituents feed back to the ecosystem and biogeochemistry. The range of values for these experiments is designed to cover and go beyond the natural range of the absorption and scattering by the water constituents. We additionally explore how different assumptions and parameterizations for the optical constituents affect the simulation results.

Detrital matter
Observations have determined that detrital matter does play a role in light attenuation, though with varying regional importance (e.g. Jerlov, 1953;Bricaud et al., 2010). We conduct several sensitivity studies to explore the relative importance of a det and b det (Fig. 17) globally in the model. We run each experiment from the same initial conditions as the "default" (EXP0) discussed in Sect. 3, and present results for the final year after 10 years of integration. We artificially alter a part det (λ) or b part det (λ) as noted below, such that a det and b det are manipulated. The experiments include the feedbacks to nutrients and productivity. In experiment EXP-D5 we explore a different parameterization for a det (λ) that was used in Fujii et al. (2007).
Removing the detrital absorption (EXP-D1) leads to bluer wavebands reaching to greater depth (Fig. 17a). This favours phytoplankton, at least in the subtropics, which absorb more efficiently in the blue part of the spectrum (i.e. Prochlorococcus, Fig. 17c) as anticipated from laboratory studies (e.g. Moore et al., 1995). On the other hand, having stronger detrital absorption (EXP-D2) leads to shallower 1 % light levels for the blue wavebands. The corresponding red-shifted light favours Synechococcus, which absorbs more efficiently in this part of the spectrum. With less irradiance absorbed in EXP-D1, we find a higher percentage is reflected at the shorter wavebands (Fig. 17d). Similarly as more irradiance is absorbed (EXP-D2), there is a reduction in the reflectance.
We observe distinct biogeochemical feedbacks. With lower absorption by detritus (EXP-D1) the depth integrated phytoplankton biomass in the high latitudes increases (Fig. 17b), leading to higher nutrient utilization in these locations. Thus the transport of nutrients to the lower latitudes is reduced (see e.g. Sarmiento et al., 2004;Dutkiewicz et al., 2005), reducing biomass in those locations. This will even further increase the 1 % light depth for the blue wavebands and consequently favour Prochlorococcus more. The lower absorption by detritus therefore leads to expansion of the oligotrophic subtropical gyres. Conversely, with more absorption (EXP-D2), we find lower depth-integrated productivity in the high latitudes, higher nutrient supply to subtropics, reduced oligotrophic regions, and favouring of Synechococcus. This feedback between the light field and the biogeochemistry can only be captured by a fully threedimensional coupled ecosystem-radiative transfer model.
The main attenuation of light with depth is through absorption, and as such alterations to the backscattering by detrital matter (EXP-D3 and EXP-D4) have little effect on the irradiance fields at depth (Fig. 17a) and thus there is little change in the dominant functional type (Fig. 17c). However scattering has a major impact on the amount and quality of the upwelling light, and as such the changes to the irradiance reflectance are large (Fig. 17d).
In EXP0, a det is calculated relative to number of detrital particles, whereas in EXP-D5 we parameterized it relative to particulate organic carbon (POC) concentrations (following Fujii et al., 2007). We find very similar patterns and magnitudes of a det (450) using these two methods. A slight difference in magnitude can be attributed to the values chosen for a POC det and p part in the respective parameterizations. There is consequently little difference in biomass, phytoplankton distributions, and reflectance between the two experiments.

Coloured dissolved organic matter
CDOM and its contribution to light absorption are observed to vary in different regions of the ocean (e.g. Jerlov, 1953;Bricaud, 1981;Nelson andSeigel, 2013, Morel et al., 2010), and many studies have attempted to empirically link a cdom to other more easily measured quantities such as Chl a (e.g. Morel, 2009). However these studies are still regional or include only sparse data. We conduct a series of sensitivity experiments that test assumptions and importance of a cdom globally and its feedback to the biogeochemistry. In two experiments (EXP-C1) and (EXP-C2) we assume no and significantly more absorption by CDOM, respectively. In additional sensitivity experiments (EXP-C3, EXP-C4, and EXP-C6) we explore the consequences of different parameterization of a cdom as used in previous model studies (e.g. Greg and Casey, 2009;Mouw et al., 2012;Fujii et al., 2007;Hickman et al., 2010).
In all experiments, a cdom (λ) is an exponential function with wavelength: In the series of experiments we make different assumptions on χ cdom : This is our default experiment detailed in previous sections.
3. EXP-C2: χ cdom = 4 · c cdom (λ o )CDOM This experiment is the same as the default (EXP0) but with CDOM artificially able to absorb 4 times as much light in each waveband.
4. EXP-C3: χ cdom = c chl (a w (λ o ) + j a chl phyj (λ o )Chl j ) Studies (e.g. Morel, 2009) have noted an empirical relationship between mean Chl a and a cdom . But regionally there is a large variation in the ratio of Chl a and a cdom (e.g. Bricaud et al., 1981;Kitidis et al., 2006;Morel et al., 2010). Here, as is done in Gregg and Casey (2007), we assume that a cdom is a function of Chl a, and c chl = 0.8 (unitless) to match the magnitudes of EXP0.

EXP-C4:
Since CDOM is part of the DOM pool, a previous model-based study (Mouw et al., 2012) assumed that some portion of the DOM pool (f cdom ) is CDOM.
Community structure shifts significantly in response to the amount of irradiance that the CDOM absorbs (Fig. 18c). No CDOM absorption (EXP0-C1) favours bluer-adapted Prochlorococcus and high absorption (EXP0-C2) leads to more Synechococcus. There is also similar impact on the biogeochemistry and shifting boundaries of the oligotrophic subtropical gyres as in the detrital experiments (Fig. 18b). The model experiments thus reveal a potentially important role for CDOM in setting phytoplankton community structure via alteration of the visible light spectrum, building on previous studies (e.g. Arrigo and Brown, 1996). The amount of absorption by CDOM impacts the reflectance, again similar to the results seen with detrital absorption (Fig. 18d).
The three alternative parameterizations of χ cdom (EXP-C3, EXP-C4, and EXP-C5) lead to very different a cdom fields (Fig. 18a). There are consequently shifts in the light fields and penetration depths of different wavebands, and corresponding regional shifts in the dominant functional type. In the parameterizations that tie χ cdom to either Chl a (EXP-C3) or DOM (EXP-C4), a cdom is almost non-existent below the 1 % light level (Fig. 18), at odds with observations (e.g. Simeon et al., 2003;Bricaud et al., 2010). Above the 1 % light level the patterns of a cdom are relatively realistic in these experiments, with higher a cdom in productive regions and lower in less productive regions. However, there are significant differences to the default run and dominant functional types are altered (Fig. 18c). The uniform a cdom simulation (EXP-C5) has a more uniform 1 % light depth along the transect, reflecting the importance of CDOM for spatial variability in the depth of the euphotic zone. Since alterations to a cdom significantly affect the irradiance propagation and hence changes the upwelling light, the impact of CDOM on the reflectance is important and all experiments show a strong response (Fig. 18d).
These experiments illustrate that the parameterization of CDOM has a very significant impact on community structure and reflectance and suggest that it is crucial to explicitly resolve CDOM in models and learn more about its variability in the ocean (Morel et al., 2010;Nelson and Siegel, 2013).

Phytoplankton
Idealized experiments were also conducted to explore the sensitivity due to phytoplankton absorption and scattering (Fig. 19). We artificially manipulate a chl phyj (λ) and b C phyj , affecting a phy and b phy .
1. EXP0: this is the default run, with each phytoplankton type having a specific absorption and scattering spectrum (Fig. 1d, e, f).
2. EXP-P1: we artificially set a chl phyj (λ) = 0 for irradiance attenuation process, but still assume that phytoplankton growth depends on light as in EXP0. This is a highly hypothetical experiment.
3. EXP-P2: we artificially set a chl phyj (λ) to 4 times that of EXP0 for irradiance attenuation process, but still assume that phytoplankton growth depends on light as in EXP0. This is therefore also a highly hypothetical experiment.
5. EXP-P4: we assume all phytoplankton have the same absorption properties (the mean, black lines, in Fig. 1d, e) for both a chl phyj (λ) and a chl psj (λ).
6. EXP-P5: we assume all phytoplankton types have the same scattering and backscattering properties (the mean, black line, in Fig. 1f).
Altering the absorption by phytoplankton (EXP-P1 and EXP-P2) has a similar impact to altering CDOM or detritus (Fig. 19). There are similar changes to the irradiance field, dominant functional type, and reflectance with consequent feedbacks to the biogeochemistry.
As discussed above, the main attenuation of light is through absorption, and thus when we assume no scattering by phytoplankton (EXP-P3) there is almost no change in dominant functional type. However, since scattering does substantially affect the upwelling light, there is some (though small) change in reflectance compared to the default run (EXP0). An experiment with 4 times b phy has similar results (not shown here).
In EXP-P4 and EXP-P5 we explore the importance of the phytoplankton type-specific absorption and scattering spectra in setting their biogeography and biogeochemical consequences. Total a phy , the irradiance field, and light penetration depths of each waveband are altered when we assume a mean absorption for all phytoplankton (EXP-P4). Total a phy is generally increased in the high latitudes and decreased at low latitudes (Fig. 19a). This occurs because diatoms (which dominate the high latitudes) have lower absorption per unit Chl a than the mean spectrum (see Fig. 1e), and picophytoplankton (which dominate the lower latitudes) have a higher absorption than the mean. Community structure is also altered ( Fig. 19c), showing that the photosynthetic absorption specific to each type is important for the emergent biogeography as has been suggested by previous studies (Bidigare et al., 1990a;Huisman and Weissing, 1995;Moore et al., 1995;Stomp et al., 2004;Hickman et al., 2010). In this study, coccolithophores have a spectrum that absorbs well in the blue-green light (Fig. 1a.) Once this advantage is removed diatoms, take over their domain. Changes to irradiance reflectance also occur as a direct result (Fig. 19d). When assuming a mean scattering spectrum for all phytoplankton (EXP-P5) we find, similar to EXP-P3, almost no difference to the irradiance field, dominant functional type, or biogeography. There are, however, small changes to the reflectance. Changes in the reflectance are also apparent when the mean a phy was used (EXP-P4). Differences in reflectance caused by phytoplankton optical properties underpin many efforts to map phytoplankton functional groups from space (see e.g. IOCCG report 15, 2014).

Discussion
In this paper we have presented a version of the MIT biogeochemistry-ecosystem model (the "Darwin Project" model) which now incorporates radiative transfer, spectrally resolved irradiance, and explicit representation of optically important water constituents. Our treatment of optical properties combines many features from prior studies (e.g. Gregg et al., 2007;Fujii et al., 2007;Mobley, 2011;Bissett et al., 1999Bissett et al., , 2004 but is more comprehensive than most. In particular, we include a detailed absorption by several different types of phytoplankton as in Gregg and Casey (2007), explicitly resolve a CDOM-like tracer as in Xiu and Chai (2014) and Bisset et al. (1999), and also resolve detrital particulate matter in a similar manner to Fujii et al. (2007).
We have evaluated our model against a range of in situ observations and satellite-derived products. The model captures the large-scale biogeochemical, ecosystem, and optical characteristics as suggested by these data sets. In particular, we have used a unique data set collected during AMT-15 which includes concurrent optical, biogeochemical, and ecosystem measurements. The model captures the observed basin scale and vertical distribution. In many of the instances where the model does not compare well to the observations, we find that the physics of the model are at least partly responsible.
The model captures spatial light absorption by different optical constituents and the relative magnitude of the scattering. However, the scattering, particularly by detrital particles, remains the least well constrained aspect (see Sect. 2.6). Each of the optical constituents resolved in the model (water, CDOM, detrital particles, and phytoplankton) has an important role in attenuating irradiance through the water column, but the relative importance differs between regions, with depth, and with wavelength (Fig. 15). CDOM was relatively more important to light absorption in highly productive regions, phytoplankton were important at the subsurface Chl a maximum and absorption by water was most important in the clear oligotrophic waters.
Our sensitivity experiments suggest that models that neglect the explicit and independently varying absorption by detrital particulate matter and CDOM are missing important components that have implications for the biogeochemistry and productivity of the model. For instance, we find that the magnitude of the light absorption of any of the water constituents that we resolve is important in setting the penetration of irradiance in different wavebands. The subsurface Chl a maximum can indeed be captured without including all constituents and spectral light (as seen in EXP-V0, and in other models; e.g. Fennel and Boss, 2003;Wang et al., 2009). However, the model developments presented were necessary for capturing the regional variability in depth of the subsurface Chl a maximum, in particular by resolving the deep pen-etration of blue-green wavelengths in the subtropical gyres. Not including any of the constituents leads to an unrealistically regionally uniform depth of the subsurface Chl a maximum.
Changes to the irradiance spectrum will have important ramifications for the community structure. Lower absorption by the optical constituents leads to deeper penetration of blue light and favours phytoplankton which absorb better in the shorter wavelengths (e.g. Prochlorococcus). However, the penetration of light also has a large impact on the biogeochemistry and biogeography at global scales. In the sensitivity studies with less light absorption, there was more primary production at the higher latitudes and reduced nutrient transport to the lower latitudes. Thus changes in absorption could impact the size of the oligotrophic regions, which in turn impact the community structure.
An important product of the model is the surface irradiance reflectance, which provides a more direct comparison to satellite data than derived products such as Chl a or primary production. These derived products rely on empirical algorithms to convert from more direct measurement of ocean colour (e.g. reflectance), which introduce a large degree of uncertainty to the output (see e.g. Campbell et al., 2002;Carr et al., 2006). Thus directly relating model output to satellite reflectance shows great promise.
The absorption by any of the optical constituents strongly determines the amount of upwelling irradiance and consequently the surface reflectance. In particular, we found that the regional variations in CDOM are important in setting the patterns of reflectance (see EXP-C5). Though alterations to scattering appears to have little effect on the in-water optical fields, they have a significant impact on the surface reflectance fields. Even slight changes to the scattering by phytoplankton (see EXP-P5) have an effect on the reflectance. Such changes are important when attempting to retrieve information on the community structure from ocean colour satellite products (e.g. IOCCG report 15, 2014).

Conclusions
The amount and type of irradiance that penetrates through the water column is an important issue when studying phytoplankton productivity and community structure. And yet, ocean models routinely offer very crude parameterizations of light attenuation and neglect the spectral quality. We have improved the MITgcm ecosystem and biogeochemistry model by incorporating spectral light, explicit radiative transfer, and representations of several optical constituents. The model performed well when compared to observations. Capturing each of the optically important constituents explicitly and including a spectrum of light was important for obtaining realistic variability in depth of the subsurface Chl a maximum and in resolving the deep penetration of blue-green wave- The sensitivity studies were intentionally hypothetical to provide a wide range of responses. They provide evidence that capturing how each of the optical constituents absorbs and scatters irradiance has important ramifications for biogeochemistry and the phytoplankton community structure. This feedback between the light field and the biogeochemistry can only be captured by a fully three-dimensional coupled ecosystem-radiative transfer model.
The model provides a platform to explore the relative importance of different optical constituents for biogeography, biogeochemistry, and optical properties such as those measured by satellite. We believe that this model will be useful in examining the role of the irradiance spectrum and pigments in setting biogeography, in exploring how changes in irradiance and/or optical constituents will impact the future oceans, and in providing a laboratory to explore the use of water-leaving radiance as a marker of changes in the marine ecosystem.

Appendix A: Ecosystem and biogeochemical model equations
The model equations are based on those of Follows et al. (2007), Dutkiewicz et al. (2009Dutkiewicz et al. ( , 2012, and Hickman et al. (2010). We consider the cycling of phosphorus, nitrogen, silica, and iron, as well as carbon, alkalinity, and dissolved oxygen (the latter three following Ullman et al., 2009). We also resolve here explicit dynamic Chl a (following Geider et al., 1998) and a tracer that mimics coloured dissolved organic matter (CDOM). We provide a complete set of the equations here.
Several nutrients, N i , nourish many phytoplankton types, P j , which are grazed by several zooplankton types, Z k . Mortality of and excretion from plankton, as well as sloppy feeding by zooplankton, contribute to a dissolved organic matter, DOM i , pool and a sinking particulate organic matter pool, POM i . Subscript i refers to a nutrient/element, j to a specific phytoplankton type, and k to a zooplankton type. Here i is PO 4 , inorganic fixed nitrogen (which includes NO 3 , NO 2 , NH 4 ), Fe, Si, and C. Particulate inorganic carbon (PIC), alkalinity (A), and dissolved oxygen (O 2 ) are also included in this framework. All tracers, X, are advected and diffused by the three-dimensional flow fields: where u = (u, v, w), velocity in physical model, K are the mixing coefficients used in physical model, and S X are sources and sinks of tracer X.
The source and sinks of each tracer, S X , are different and including biological transformations, chemical reactions, and external sources and sinks. Phytoplankton are assumed to have fixed elemental ratios following Redfield (1934). The base currency of the plankton equations is phosphorus. Nutrients: Plankton: Chlorophyll a: Particulate and dissolved matter: Alkalinity: Dissolved oxygen: where µ j is the growth rate of phytoplankton j (function provided below), M ij is the matrix of ratios of element i to phosphorus for phytoplankton j , r dom i is remineralization rate of DOM for element i, here P, Fe, N, C, r pom i is degradation/remineralization rate of POM for element i, here P, Si, Fe, N, C, d cdom is degradation rate of CDOM to DOM for element i, here P, Fe, N, C, γ T is temperature regulation of biological rates (function provided below), c scav is scavenging rate for free iron (function provided below), Fe is free iron (description provided below), F atmos is atmospheric deposition of iron dust on surface of model ocean, F sed is the sedimentary source of iron (function provided below), ζ no3 is oxidation rate of NO 2 to NO 3 (function provided below), ζ no2 is oxidation rate of NH 4 to NO 2 (function provided below), no3 j is fraction inorganic nitrogen uptake from nitrate (function provided below), no2 j is fraction inorganic nitrogen uptake from nitrite (function provided below), nh4 j is fraction inorganic nitrogen uptake from ammonium (function provided below), H ocrit = 1 if O > O crit and 0 if O =< O crit , O crit is critical oxygen level for denitrification, R denit is N : P ratio in denitrification, R dno3 is ratio of NO 3 relative to all N in denitrification, D denit is denitrification rate (function provided below), R rj is ratio of inorganic carbon to organic phosphorus produced by phytoplankton j , F C is air-sea flux of carbon dioxide (function provided below), D C is dilution/concentration of carbon by addition/loss freshwater, D A is dilution/concentration of alkalinity by addition/loss freshwater, F O 2 is air-sea flux of oxygen (function provided below), d pic is dissolution rate of PIC, m pj is mortality/excretion rate for phytoplankton j , m zk is mortality/excretion rate for zooplankton k, m z2k is quadratic mortality for zooplankton k, g j k is grazing of zooplankton k on phytoplankton j (function provided below), ζ j k is grazing efficiency of zooplankton k on phytoplankton j (function provided below), w pj is sinking rate for phytoplankton j , w pom i is sinking rate for POM i, w pic is sinking rate for PIC, ρ j is Chl a : C of new growth (function provided below), θ j is local Chl a : C ratio, θ oj is acclimated Chl a : C (function provided below), t chl is acclimation timescale for Chl a, ϕ mp ij is fraction of dead/respired phytoplankton organic matter that goes to DOM i , ϕ mz ik is fraction of dead/respired zooplankton organic matter that goes to DOM i , ϕ g ij k is fraction of sloppy grazing that goes to DOM i , f cdom is fraction of DOM produced that enters CDOM pool, ι cdom is bleaching rate for CDOM, λ=700 λ=400 E 0 (λ) is local total scalar irradiance, I cdom is PAR above which CDOM bleaches.

A1 Temperature regulation of biological rates
Biological rates (plankton growth and the parameterization of remineralization of organic matter) are represented as a function of temperature, following the Arrhenius equation (Kooijman, 2000), similar to Eppley (1972): where τ 1 is a coefficient to normalize the maximum value, A E and T o regulate the form of the temperature modification function, and T is the local model ocean temperature.

A2 Phytoplankton growth
Phytoplankton growth is a function of temperature, irradiance, and nutrients. We follow Hickman et al. (2010), which in turn follows Geider et al. (1998), such that the growth rate is equal to the carbon-specific photosynthesis rate: where P C mj is light-saturated photosynthesis rate (see function below), Ej is the scalar irradiance absorbed by each phytoplankton multiplied by φ max j , the maximum quantum yield of carbon fixation (see function below), and θ j is Chl a : C for each phytoplankton (see function below).
The light-saturated photosynthesis rate is a function of nutrients and temperatures: where P C m max j is maximum photosynthesis rate of phytoplankton j , γ T is modification of growth rate by temperature (see above), and γ Nj is modification of growth rate by nutrients for phytoplankton j (see function below).
For each phytoplankton j , is the scalar irradiance absorbed by each phytoplankton, j , multiplied by φ max j , the maximum quantum yield of carbon fixation (units mmol C (mg Chl a) −1 d −1 ). This is by definition the spectrally resolved product of irradiance and the initial slope of the Chl a normalized photosynthesis versus irradiance curve. a chl psj (λ) is the Chl a-specific photosynthetic absorption spectra in each waveband λ.
The local Chl a : C ratio θ j is The increase of Chl a due to growth term (M Cj ρ j µ j P j ) in Eq. (A11) follows Geider et al. (1998), with and the acclimated Chl a : C follows Geider et al. (1997): where θ maxj is the maximum Chl a : C ratio each phytoplankton can reach. Phytoplankton can be photoinhibited (following Hickman et al., 2010), such that P C j reduces to P C inhib j above E kj : where κ inhib is the inhibition coefficient and E kj is the light saturation parameter.
where a chl psj (λ) is the mean light absorption by photosynthetic pigments between 400 and 700 nm. Nutrient limitation is determined by the most limiting nutrient: Limitations by PO 4 , Si, and Fe are all parameterized following the Michaelis-Menton formulation: where κ N ij is the half-saturation constant of nutrient i = PO 4 , Si, Fe for phytoplankton j . Nitrogen is available in three forms, of which ammonia is the preferred type: where κ in j is the half-saturation constant of IN = NO 3 + NO 2 , κ nh4 j is the half-saturation constant of NH 4 , and ψ reflects the fixed nitrogen uptake inhibition by ammonia.

A3 Zooplankton parameterization
Zooplankton grazing is parameterized as where g max j k is maximum grazing rate of zooplankton k on phytoplankton j ; η j k is palatability of plankton j to zooplankton k; G k is palatability (for zooplankton k) weighted total phytoplankton concentration, equal to j [η j k P j ]; κ pk is the half-saturation constant for grazing of zooplankton k; and n is exponent for Holling type II or III (n = 1 or 2), in this study n = 2. The maximum grazing g max j k depends of the relative size of the phytoplankton j and zooplankton k, with a faster rate if they are both small or both big (g max a ), and slower if they are in different size classes (g max b ).
Zooplankton are assumed to have both a linear and quadratic loss term. The linear term represents excretion and mortality; the quadratic loss terms represent grazing by higher trophic levels (Steele and Henderson, 1992) that are not explicitly resolved in this model.

A4 Nitrogen cycle
Phytoplankton take up DIN in three forms (NH 4 ,NO 2 ,and NO 3 ). To separate out how much comes from each source we have the functions in Eqs. (A5)-(A7): The oxidation of NH 4 to NO 2 and NO 2 to NO 3 is parameterized as a function of the total scalar irradiance: where ζ ono3 and ζ ono2 are maximum rates and I 0 is critical light level below which oxidation occurs. Denitrification occurs when O < O crit , in which case O 2 is not used during remineralization but instead NO 3 is used, such that We assume the denitrification formula suggested by Anderson (1995) for determining R denit : C 117 N 16 P + 120NO 3 ⇒ 117CO 2 + PO 4 + 68N 2 .

A5 Iron parameterization
The iron model we use is based on that of Parekh et al. (2004Parekh et al. ( , 2005. We explicitly model the complexation of iron with an organic ligand: where Fe and L are free iron and ligand, respectively, FeL is ligand-bound iron, L T is total organic ligand (assumed to be a constant), β fe = k f k d is ligand binding strength, k f is the forward rate constant, and k d is the reverse rate constant.
We assume that only the free iron (Fe ) can be scavenged, c scav Fe , and parameterize this as a function of the particulate organic carbon (POC) present (empirical values based on those found for thorium; Honeyman et al., 1988); a similar approach was used in Parekh et al. (2005): where c o determines maximum scavenging rate for iron, ξ is an empirically determined constant, and R C : P is the carbon to phosphorus ratio of the POM. The sedimentary source (F sed ) is parameterized as a function of the sinking organic matter reaching the ocean bottom as suggested by Elrod et al. (2004): where R sed is the ratio of sediment iron to sinking organic matter.

A6 Air-sea exchange
Air-sea exchange of CO 2 and O 2 are given by [CO 2 ] is calculated from DIC and alkalinity concentrations following Follows et al. (2006), which included deducing the pH at all surface locations. The gas transfer coefficient is parameterized following Wanninkhof (1992) and is a function of the wind speed and Schmidt number (a function of sea surface temperature). [CO 2 ] sat is determined as a function of partial pressures of CO 2 in the air, atmospheric pressure, sea surface temperature, and salinity. [O 2 ] sat is provided by Garcia and Gordon (1992). All coefficients of the air-sea flux calculations are determined using the algorithms used in the Ocean Carbon-Cycle Model Intercomparison Project (OCMIP) (e.g. Matsumoto et al., 2004). The radiance in the ocean in its most general form, L(x, θ, ϕ, λ), depends on location and orientation in addition to wavelength (units W m −2 sr −1 nm −1 ). Neglecting horizontal gradients, the z dependence of L is described by the classical radiative transfer equation, dL(θ, ϕ) dz cos θ = −cL(θ, ϕ) + β(θ, ϕ, θ , ϕ )L(θ , ϕ )d , where β(θ, ϕ, θ , ϕ ) is the rate of scattering of light from θ , ϕ into θ, ϕ. We assume that the ocean is optically isotropic, so β is invariant under simultaneous rotation of original and scattered angles (in fact it depends only on the relative angle). The integral over one set of angles therefore yields an angle-independent value, 4π β(θ, ϕ, θ , ϕ ) d = 4π β(θ, ϕ, θ , ϕ ) d = b.
Here, b is then the total scattering coefficient and the total scattered light is β(θ, ϕ, θ , ϕ )L(θ , ϕ ) d d = b L(θ , ϕ ) d = bE 0 and may be decomposed into forward and backward scatter- The attenuation coefficient c represents loss due to absorption and scattering, c = a + b.
At the sea surface, the downward part of L(θ, ϕ) for θ < π/2 is required to equal the output of the atmospheric radiative transfer model (OASIM). The ocean is assumed to be infinitely deep, with vanishing light at infinite depth.
In terms of the effective backscattering coefficients, In general,ῡ s andῡ u depend on the angular profile of the radiation field, and r s and r u , which describe the scattering of downward into upward and upward into downward radiation, depend on both the scattering function and the radiation field. We close the system of equations by making the following assumptions (following Aas, 1987): Equations (B3)-(B5) are the three-stream equations (given in main text as Eqs. 1-3, though note that here we dispense with function of λ for simplicity).
The equation for E d (Eqs. B3 or 1) is readily integrated, In contrast to Aas (1987), Ackelson et al. (1994) and Gregg (2002) we do not make further approximations but instead solve the remaining equations explicitly. We can write the remaining two equations (Eqs. B4 and B5, as well as Eqs. 2 and 3) as where and M, F d , and B d are assumed to be piece-wise constant as a function of z.

S. Dutkiewicz et al.: Modelling optical properties
Following Kylling (1995) we write the inhomogeneous solution as where x and y satisfy the equation with solution The eigenvalues of M are Within a computational layer, the general solution can be written as where r + = R 2 = B s /D and r − = 1/R 1 = B u /D. The offsets in the exponents have been introduced so that both exponentials are smaller than 1. The coefficients c + and c − have to be determined from boundary conditions. At the sea surface, we require that E s and E d coincide with the output of OASIM, In the bottom layer, k bot , we require zero light at infinite depth, i.e. c − k bot = 0. At layer boundaries, z k+1 , we require continuity, In order to solve this coupled system of equations, we follow Kylling et al. (1995) and Toon et al. (1989), who observed that it can be transformed to tri-diagonal form by eliminating c − k+1 from the first equation, and c + k from the second equation, where e + k = e −κ + k (z k+1 −z k ) and e − k = e −κ − k (z k+1 −z k ) . The reduced system is solved explicitly using Gaussian elimination. In order to estimate backscattering b b from the observations made during AMT-15 we utilize the measured downwelling irradiance, E dn , and upwelling, zenithward radiance, L u . We use the procedure of Boynton (1997, 1998) with the radiative transfer package DISORT, version 2.0 beta. We use the Boynton (1997, 1998) parameterization rather than the quasi-analytical algorithm (Lee et al., 2002(Lee et al., , 2007 since we are dealing with profiles and not surface water-leaving radiance. Gordon and Boynton (1998) propose that R = E u /E dn and X = b b /a are related as We drop E dn (z) from numerator and denominator and discretize as 2k j and k j = 1 z j +1 − z j ln E j E j +1 .
In order to solve for X, we write and get For noisy data, this estimate of X may become negative. We drop the derivative term where this happens -i.e. X is approximated by 3R.