Optical Modeling of Spectral Backscattering and Remote Sensing Reflectance From Emiliania huxleyi Blooms

In this study we develop an analytical model for spectral backscattering and ocean colour remote sensing of blooms of the calcifying phytoplankton species Emiliania huxleyi. Blooms of this coccolithophore species are ubiquitous and particularly intense in temperate and subpolar ocean waters. We first present significant improvements to our previous analytical light backscattering model for E. huxleyi coccoliths and coccospheres by accounting for the elliptical shape of coccoliths and the multi-layered coccosphere architecture observed on detailed imagery of E. huxleyi liths and coccospheres. Our new model also includes a size distribution function that closely matches measured E. huxleyi size distributions. The model for spectral backscattering is then implemented in an analytical radiative transfer model to evaluate the variability of spectral remote sensing reflectance with respect to changes in the size distribution of the coccoliths and during a hypothetical E. huxleyi bloom decay event in which coccospheres shed their liths. Our modeled remote sensing reflectance spectra reproduced well the bright milky turquoise colouring of the open ocean typically associated with the final stages of E. huxleyi blooms, with peak reflectance at a wavelength of 0.49 μm. Our results also show that the magnitude of backscattering from coccoliths when attached to or freed from the coccosphere does not differ much, contrary to what is commonly assumed, and that the spectral shape of backscattering is mainly controlled by the size and morphology of the coccoliths, suggesting that they may be estimated from spectral backscattering.


INTRODUCTION
Coccolithophores are phytoplankton that form an exoskeleton of calcite scales called coccoliths (Figure 1). They are major oceanic calcite producers and carbon exporters in the ocean (Iglesias-Rodríguez et al., 2002;Broecker and Clark, 2009). The cosmopolitan coccolithophore species Emiliania huxleyi (Lohmann) Hay & Mohler thrives from tropical to subpolar oceans and forms large-scale, intense blooms particularly in temperate waters (Tyrrell and Merico, 2004). In the later bloom stages E. huxleyi overproduces and sheds coccoliths, giving the ocean a bright milky-turquoise appearance, which makes these blooms easily discernible from optical satellites ( Groom and Holligan, 1987;Holligan et al., 1993;Brown and Yoder, 1994), especially ocean color satellites Balch et al., 2005). The backscattering of light by coccoliths and coccospheres forms the basis of the remote sensing algorithm for the retrieval of the concentration of particulate inorganic carbon, PIC (see Table 1 for a list of symbols and abbreviations), or calcite Balch et al., 2005). Remotely sensed PIC has been widely used to study coccolithophore bloom extent, occurrence, and timing in the global ocean (Balch et al., 2011;Hopkins et al., 2015), as well as their poleward expansion in response to climate change (Winter et al., 2014;. In situ and laboratory observations show that the relationship between light backscattering and PIC may vary, however, and the causes for this variability are poorly understood (e.g., Balch et al., 1991). While it is commonly assumed that PIC-specific light backscattering (=backscattering per unit PIC) is orders of magnitude higher for detached coccoliths than for coccospheres (Balch et al., 1991(Balch et al., , 1996, Gordon et al. (2009) suggested that this would not be so based on a comparison of an optical model and in situ measurements in a coccolithophore bloom. Optical models are indeed useful to investigate variations in light backscattering due to changes in size and morphology of coccoliths and coccospheres, but investigations have been limited due to the computational demand of commonly used optical models suitable for light backscattering such as the discrete dipole approximation (Gordon and Du, 2001;Gordon, 2006;Zhai et al., 2013). Recently, however, Fournier and Neukermans (2017) developed an analytical optical model for light backscattering from coccoliths and coccolithophores of E. huxleyi. Our previous model closely matched results obtained from the discrete dipole approximation computations of Zhai et al. (2013) in which E. huxleyi coccoliths were represented as two thin circular disks attached together by a tube and perforated by a set of wedge shaped openings.
In the present paper we first improve our previous analytical optical model to account for a more realistic representation of the morphology of E. huxleyi coccoliths and coccospheres. More specifically we account for (i) the fact that coccoliths are elliptical in shape with triangularly shaped wedge openings (Figure 1), and (ii) detailed coccosphere architecture with interlocking coccoliths as obtained from three dimensional electron microscopy observations of E. huxleyi coccolithophores (Hoffmann et al., 2015). We then improve our previous size distribution model based on detailed analyses of observed coccolith size distributions (Young et al., 2014). Using our new coccolithophore and coccolith backscattering model and an algebraic radiative transfer model we investigate the variability in spectral backscattering and ocean color remote sensing reflectance (i) due to changes in E. huxleyi coccolith size distribution and morphology, (ii) due to changes in the ratio of free to attached coccoliths, and (iii) during a hypothetical bloom decay event in which E. huxleyi sheds its attached liths. Programming codes for all computations in this paper are provided in Neukermans and Fournier (2018), including future releases.

METHODS
The light backscattering coefficient of calcite particles suspended in seawater, b bpic (in m −1 , see Table 1 for a list of symbols and units), can be partitioned in its contributions from liths and coccospheres, represented by b bl and b bc , respectively: For populations of identical coccoliths and coccospheres we can write: b bpic = N l σ bbl + N cc σ bbs where N l and N cc are the number concentrations of coccoliths and coccospheres (in m −3 ), respectively, and σ bbl and σ bbs are the corresponding backscattering cross-sections for individual a o , a r , a w , a t Components of a, µm (see Figure 2) a om Minimum semi-major axis length for a distribution of coccoliths, µm a peak Location of the peak value of the shifted gamma particle size distribution, µm Mean backscattering cross-section for a population of elliptical coccoliths, m 2 σ g Geometric cross-section of a particle, m 2 σ gl , σ gs Geometric cross-section of a coccolith or coccosphere, m 2 σ gl Geometric cross-section averaged over the coccolith size distribution, m 2 τ l Thickness of the elliptical coccolith, µm τ n Dimensionless thickness parameter of the elliptical coccolith w r , w w , w t Width of rings comprising the elliptical coccolith, µm (see Figure 2) w c Critical width of the ring where the gap between the pillars of the coccolith becomes equal to λ w /4, µm x (λ) Backscattering albedo, dimensionless ξ Dimensionless spectral slope parameter of b bp (λ) coccoliths and coccospheres (in m 2 ). The backscattering crosssection is defined as the product of the geometric cross-section of the calcite particle σ g (in m 2 ) and the dimensionless efficiency factor for backscattering of the particle Q bb so that: σ bbl = Q bbl σ gl and σ bbs = Q bbs σ gs where the efficiency factor Q bb is the ratio of radiant power backscattered by the particle to radiant power intercepted by the geometric cross-section of the particle (e.g., Morel and Bricaud, 1986). The underlying assumptions to Equation (3) are, first, that coccoliths and coccolithophores are randomly oriented with respect to the direction of incident light. This allows us to evaluate σ g through the theorem of Cauchy (1850) which states that the average σ g of a randomly oriented particle of convex shape is equal to a quarter of its surface area. Second, for coccolithophores and coccoliths made of high refractive index calcite material, the backscattering efficiency Q bb is dominated by the contribution of reflection from the particle surface (Fournier and Neukermans, 2017). For an ensemble of randomly oriented reflecting surfaces the theorem of van de Hulst (1957) formally states that: "the [angular] scattering pattern caused by reflection on large convex particles with random orientation is identical with the [angular] scattering pattern of a large sphere made of the same material and with the same surface condition." This equivalence is a direct consequence of the angular averaging procedure, which holds for randomly oriented surfaces. Therefore, it follows that the unity-normalized angular distribution of reflection from a randomly oriented convex body is shape-independent and can be characterized by a single scattering phase function in the back-direction. This allows us to treat orientation averages of the scattering efficiency factor and the geometric cross-section separately, which is implicitly assumed in Equation (3). In short, for randomly oriented highly reflective particles such as coccoliths and coccolithophores, the average σ bb can be obtained by multiplying σ g with Q bb .
In what follows we first develop expressions for σ bbl for single coccoliths (section Model for the Backscattering Cross-Section of an Elliptical Coccolith) as a function of lith size, shape, and morphology. We thereby build on our previous modeling work (Fournier and Neukermans, 2017) which relies on the fundamental hypothesis that for particles large compared to the wavelength of light, the backscattering is dominated by the reflection from the surfaces of the particles. This approximation holds particularly well for coccoliths because they are made of calcite, a material with a refractive index of 1.2 relative to water and therefore a high reflection coefficient. The reflection from the four calcite surfaces which comprise the coccolith particle is considered to be specular, and the surfaces are assumed smooth when the surface roughness features are smaller than a quarter of the light wavelength in water, λ w . If the surface roughness features become larger than λ w /4, light is reflected diffusely. This is represented in our previously derived formula for the backscattering efficiency of a single circular coccolith particle: where F rg is the maximum fraction of the projected area of the particle that can become rough, RF rg is the actual fraction of the area of the particle that contributes to rough scattering, ω bsm represents the specular reflection component coming from the smooth part of the particle surface and ω brg is a diffuse reflection component originating from the rough part of the particle surface. Expressions for ω brg and ω bsm as functions of the refractive index are derived in Fournier and Neukermans (2017).

Model for the Backscattering Cross-Section of an Elliptical Coccolith
To model the light backscattering cross-section of an elliptical coccolith of E. huxleyi we need to derive Q bbl and σ gl (Equation 3) and therefore we need to consider the detailed geometry of a coccolith (Young et al., 2014) which is schematically depicted in Figure 2. The semi-major axis a is ½ the coccolith length while the semi-minor axis b is ½ the coccolith width. As indicated in Figure 2, the structure of the coccolith can be accurately represented by a set of four ellipses separated by three rings of fixed width. We have an outermost ellipse separated by a thin ring of width w r from the ellipse which marks the beginning of the zone with the open wedges. This ellipse is in turn separated by a ring of width w w from an inner ellipse. These ellipses define the boundaries of the lith area containing the wedged openings. The inner ellipse is separated from a center ellipse by a width w t that represents the boundaries of the inner tube joining the distal and proximal sheets. The average geometric cross-section of a randomly oriented elliptical coccolith is derived in Supplementary Information section 1.1 (based on Cauchy, 1850's theorem) and repeated here for convenience: where τ l is the thickness of the elliptical disk and P e a, b its perimeter with equations given in Supplementary Information section 1.1. We can represent the proportional relations between the various ellipse parameters (a o , a r , a w , a t ) and b o , b r , b w , b t as follows (Figure 2): For the coccolith in Figure 1, with major axis in µm of a o = 1.8688, we have the following proportionalities: The maximum fraction of the surface area of the coccolith that can become rough F rg lies between ellipses with parameters (a r , a w ) and b r , b w so we can write: which gives F rg = 0.5418 for the coccolith in Figure 2.
We will now analyze more closely the structure of wedges and pillars. The structure of the wedges and pillars in the zone between a r and a w can be represented as a set of gaps with maximum widths s gmax and minimum pillar width s pmin respectively along the elliptical boundary at a r and minimum widths s gmin and maximum pillar width s pmax . The gap width narrows down to s gmin on the elliptical boundary at a w while the pillar width grows to fill the empty space. Given N p pillars around the circumference of the ellipse at a r we have: N p s gmax + s pmin = P e a r , b r N p s gmin + s pmax = P e a w , b w Defining the ratio of the minimum pillar width to the maximum gap width, r pg , and the ratio of minimum to maximum gap width, r sh , as: r pg = s pmin s gmax and r sh = s gmin s gmax (9) we obtain the following expression for s gmax where P n is the dimensionless ellipse perimeter given in Equation (S13). As we proceed from the ellipse at a r to the ellipse at a w the gap s g narrows from s gmax to s gmin as w goes from 0 to w w To determine where the backscattering becomes rough we need to compute at what value of a the gap s g becomes equal to 1/4 of λ w (Gordon, 2006;Fournier and Neukermans, 2017). We define w c as the critical width where the gap becomes equal to λ w /4, then Equation (11) gives: The corresponding critical semi-major axis a c and semi-minor axis b c are given by: As long as the gap s g is smaller than s gmax and larger than λ w /4 we can compute the critical width as: When s g is smaller than λ w 4 we must set w c to 0 and on the high side as the wavelength gets smaller we must limit w c to a maximum value of w w : We can then proceed to compute the ratio w c a o which is required to evaluate the fraction of rough surface of the coccolith using Equations (15) and (10).
Substituting this expression in Equations (13) and (14) we obtain directly the fraction of the elliptical surface of the coccolith that is rough: The only variable determining the behavior of RF rg is the ratio of the wavelength in water λ w to the major axis a o of the elliptical coccolith. The rest of the elements are fixed by the proportionality factors dictated by the assumed shape of the coccolith. Combining Equations (8) and (9) we find for the ratio of circumferences of the outer and inner ellipse where the wedges and pillars occur: is the ratio of maximum to minimum pillar width. For triangularly shaped wedges such as those measured by Young et al. (2014) s gmin = 0 = r sh and using Equation (S13) for a non-dimensional version of the perimeter ratio, P n , we find: Equation (19) then simplifies to: For the coccolith morphology studied by Young et al. (2014) and shown in Figure 1 we have r pg = 0.666 and f pb = 1.602 which means that the base of the pillar s pmax is almost as wide as the wedge opening at the top s gmax . We will use these values in the present work as it gives the best fit to the structure of the coccolith pillars of the distal sheet.

Backscattering Cross-Section for a Population of Elliptical Coccoliths
To evaluate the geometric cross-section averaged over a size distribution we rewrite Equation (5) in a non-dimensional form (derived in Supplementary Information section 1.1, Equation S14): with P n the dimensionless ellipse perimeter defined in Equation (S13). We can then obtain the resulting mean backscatter cross section by multiplying with the backscatter efficiency, as in Equation (3). For convenience we define the thickness parameter as a non-dimensional proportionality: which gives In our previous work (Fournier and Neukermans, 2017) we assumed a first-order shifted gamma size distribution P (a o ) of coccoliths. Detailed analyses of the shape of measured coccolith size distributions (Young et al., 2014), however, show that these are better approximated by a second-order shifted Gamma function (Supplementary Information section 1.2) of the following form: where a om is the minimum size of the coccoliths. The following relations can be used to express the parameters in terms of the measured mean µ and standard deviation σ of the coccolith populations: where a peak is the location of the peak value of the distribution. Using Equation (25), the geometric cross-section averaged over the coccolith size distribution, σ gl , is therefore: The mean backscattering cross-section and backscattering efficiency factor for a population of elliptical coccoliths are then given by, respectively: Note that P (a o ) can now be expressed as a pure function of measured mean µ and standard deviation σ by using the relations in Equation (27): Coccolithophore Architecture To address whether light backscattering from E. huxleyi coccospheres differs in any way from the backscattering by the sum of its individual coccoliths we need to look into the details of coccolithophore architecture. In our previous work (Fournier and Neukermans, 2017) we assumed that coccospheres were composed of coccoliths which completely cover the sphere in multiple layers attached on top of the core and on top of one another. This determined the number of coccoliths attached to a given core for a given number of layers. Three dimensional electron microscopy observations of coccolithophores (Hoffmann et al., 2015), however, suggest an entirely different layering of coccoliths on coccospheres, which we expect to impact the light scattering properties of coccolithophores and the number of liths comprising the coccolithophore, and by consequence the difference in backscattering from attached vs. free liths. We will first present the coccosphere model proposed in our earlier work (Fournier and Neukermans, 2017), referred to as the stacked liths coccosphere model, and next present an improved coccosphere architecture model, referred to as the interlocking liths coccosphere model.

Stacked Liths Coccosphere Model
In our previous work we assumed that E. huxleyi coccolithophores were built up of coccoliths which completely cover the coccosphere in multiple layers with liths stacked on top of one another as shown in Figure 3A. As per Equation (3), we examine the geometric-cross sections corresponding to E. huxleyi coccoliths when attached to and freed from the coccosphere. As the average geometric cross-section of any FIGURE 3 | Schematic representation of (A) stacked and (B) interlocked liths coccosphere models. Three dimensional electron microscopy images indicate that the interlocked pattern is the one closest to reality (Hoffmann et al., 2015). randomly oriented convex shape is equal to ¼ of its surface area (Cauchy, 1850) we find for a coccosphere of radius r s that For a randomly oriented elliptical disk of semi-major axis a, semiminor axis b, aspect ratio r = a/b, and thickness τ d , the mean geometric cross-section is: where f d = τ d /a, the disk thickness relative to a. Note that we have used the first order approximation π a + b for the outer perimeter of the ellipse P e a, b as this keeps the algebra simple and is accurate enough for coccoliths who are nearly spherical (r = 1.2) (see Equation S7 in the Supplementary Notes). Based on Paasche (2001), we assume that in the first layer the coccoliths' proximal sheets are attached to the core and completely cover it, which gives the following equality: where N is the number of liths covering the sphere. Note that the simple elliptical disk model implies that the surface of the proximal and distal sheets is the same. If the core of the sphere is not absorbing, the liths on both the front and back surfaces of the sphere contribute to the backscattering, i.e., Q bbs = 2Q bbl , and the total backscattering cross-section of the coccosphere is given by: Using Equations (33) and (34), the backscattering cross-section for N disks, σ bbNl , is given by: which can be simplified to: From Equations (35) and (37), we thus find that for a nonabsorbing coccosphere covered with a single layer of N liths the ratio of the backscattering cross-section for free liths to that of the coated coccosphere is: This expression is near unity for thin disks (f d ≈ 0). To account for the multi-layered structure of the coccosphere, let O vl be the number of layers. Adding (O vl − 1) additional layers of coccoliths attached by their proximal sheets to the distal sheets of the coccoliths on the layer beneath them we have: where r s is the sphere radius increase due to each layer of liths. If the N proximal sheets also completely cover each upper layer we obtain the following general formula for the ratio of the backscattering cross-section for free liths to that of the coccosphere with O vl number of layers: We thus find that the backscattering from non-absorbing coccospheres with coccoliths attached is only slightly smaller than the backscattering from all its individual liths at any given wavelength. This near-equality is due the fact that the geometric cross-sections of randomly oriented disks covering a sphere and of the covered sphere are equal if one assumes the sphere is transparent and the reflections from both the front and back layers are the dominant contributors to the backscattering as is the case for calcite coccoliths with their high index of refraction. When the core of the coccosphere is absorbing, the area from the back surface shadowed by the core no longer contributes to the backscattering from the coated coccosphere and we have: (Fournier and Neukermans, 2017): Where Q as k, 2r c is the two-way absorption efficiency, i.e., the absorption efficiency for light going to the back surface and coming back toward the front of the core after reflection from the back surface. Using the anomalous diffraction approximation for spheres the two-way absorption efficiency of the spherical core is (Fournier and Neukermans, 2017): To evaluate Q as k, 2r c we use the spectral coccolithophore core absorption k values obtained by Frontiers in Marine Science | www.frontiersin.org Stramski et al. (2001) and Neeley et al. (2015). Using Equation (39), we can now rewrite Equation (41) So in the case of thin liths (f d ≈ 0) and layers (r s ≈ 0) and a fully absorbing core (Q as k, 2r c ≈ 1) this would imply that the scattering from the free liths would be twice as much as the scattering from the coated coccosphere: We thus find that backscattering from fully absorbing coccospheres with coccoliths attached is at least half as much as the backscattering from all its individual liths. Therefore, at wavelengths where the core absorbs most, i.e., in the chlorophylla absorption peaks at 0.44 and 0.68 µm, the ratio in Equation (44) will be higher, but will not reach 2 as E. huxleyi cells are not perfect black bodies (Morel and Bricaud, 1981).

Interlocking Liths Coccosphere Model
Three-dimensional ion beam microscopy images of E. huxleyi coccolithophores indicate that coccoliths are actually interlocking so that the central tube of the coccolith penetrates through the layer below, such as shown in Figure 3B (Hoffmann et al., 2015). The distal sheets completely cover the outer surface of every lith layer and the proximal sheets overlap thus reducing the area available and the number of coccoliths that can be accommodated per layer. A further implication is that, to first order, each layer of coccoliths contains about the same number of liths. How many liths can be accommodated on a coccosphere by the interlocked liths model and how does this compare to the amount accommodated by the stacked liths model? Accounting for the fraction of the space occupied by the coccolith tube structure which has a semi-major axis a w , we find the following equality: where a 0 is the semi-major axis of the outer part of the coccolith, r c is the radius of the coccolithophore core, and N il is the number of coccoliths on the first layer. Using Equation (34) to determine the number of liths per layer that can be accommodated by the stacked liths model, N sl , we find the ratio of the number of liths per layer for the two models: Expressing the exclusion factor f 2 t in terms of the geometry of the coccolith we have: For liths with the morphology measured by Young et al. (2014) we have Measurements on cultured E. huxleyi cells give typical values for the mean core radius r c of 2.2 µm and the mean semi-major axis a 0 of 1.25 µm (Hoffmann et al., 2015), which are at the lower end of the values measured in situ (Young et al., 2014). The ratio a 0 /r c of 0.56, however, is close to the value of 0.46 used by Zhai et al. (2013) and us (Fournier and Neukermans, 2017) in our model studies for circular coccoliths. Accounting for the elliptical shape of coccoliths using the equivalent radius a eq : a eq = a 0 b 0 = 1.25 × 1.25/1.2 ≈ 1.14 we find that a eq /r c is 0.52, which brings the coccolithophore model of Zhai et al. (2013) even closer to the experimental architecture of E. huxleyi coccospheres and liths observed by Hoffmann et al. (2015). We will therefore use the same proportions here as found in the model of Zhai et al. (2013) for the various other parameters of the liths such as distal and proximal sheet thicknesses (τ l /r c =0.07), the separation between the sheets (d h /r c =0.18), and the tube width (a t /r c =0.18). Solving Equation (34)  To determine the number of coccoliths attached on additional layers of the coccosphere, we use the measurements of diameter increment per layer obtained from laboratory cultures, which are on the order of 1.0 µm (Hoffmann et al., 2015). For the stacked liths coccosphere model this increase in diameter per layer should be larger by twice the sum of the thicknesses of the proximal and distal sheets of the coccolith. Using the same geometry as in our previous paper (Fournier and Neukermans, 2017), this implies that the diameter increase per layer, r s , should be 1.53 µm. For l layers this gives N sl = 4 r r c + l r s 2 a 2 0 (50) Table 2 compares the number of liths for multilayered coccospheres between both coccosphere models. We immediately see that the interlocked liths model gives rise to many less liths being attached to a given core for the same number of layers, which is more consistent with experimental measurements (Hoffmann et al., 2015) than what we obtain from the stacked liths model. The analysis of the ratio of the backscattering cross-section for free liths to that of the coated coccosphere with interlocked liths is similar to that for the stacked liths model. The total area of the core covered by the liths stays the same for all layers and therefore the total lith area is: We thus find the following ratio for a non-absorbing core: which is exactly the same as for the stacked liths model given in Equation (40). For an absorbing core, we also find the same ratio as for the stacked liths model: Only the absolute magnitudes of the cross-sections between both models are different since for the first layer we have: Based on the coccolithophore structure measured by Hoffmann et al. (2015), the lith morphology obtained by Young et al. (2014), and our previous results we determine the following general values for the key parameters in Equation (53): In summary, no matter which coccosphere architecture model we use, we get the fundamental result that backscattering from detached coccoliths is not so much greater than the backscattering from the coccolithophore with its coccoliths attached, which confirms the suggestion of Gordon et al. (2009), and that the difference in backscattering is governed by the absorption of the core which has a spectral signature.

Changes in Light Backscattering and Remote Sensing Reflectance During E. Huxleyi Bloom Decay
We simulate changes in spectral backscattering during a hypothetical E. huxleyi bloom decay phase in which E. huxleyi sheds all its attached liths. We will account for two possibilities for the fate of the E. huxleyi cores; they either die during bloom decay (for example due to viral lysis), or they remain without liths . Their presence only affects light absorption (and by consequence remote sensing reflectance) as we assume that the calcite material dominates the backscattering coefficient, in agreement with observations in natural E. huxleyi blooms. These bloom dynamics scenarios are based on observations of bloom decay in laboratory experiments and in situ  with the following simplifications: (i) coccolith production stops when shedding begins so that the number of liths remains constant during the process and none are detached at the beginning, and (ii) dead E. huxleyi cores do not absorb light and the light absorption properties of the cores that remain without liths remain constant throughout the bloom event. Thus, initially we have a suspension of living E. huxleyi cells, coccospheres with absorbing cores and with all N ltot liths attached. This suspension then gradually transforms into a suspension comprised of free liths with a certain fraction of surviving absorbing cores left, F sc , which varies between 0 and 1. If we define F nc as the ratio of cores that have shed their liths to the total initial number of coated cores (F nc = number of naked cores/N cc ), which varies from 0 to 1 during the coccolith shedding event, we can express the total backscattering cross-section σ bbtot at any given point during the event as follows: with the backscattering cross-section of a coated coccolithophore.

Algebraic Radiance Model
To evaluate variability in ocean color remote sensing due to changes in the size distribution of E. huxleyi coccoliths and during bloom decay, we develop an algebraic remote sensing reflectance model. The model is based on the work of Albert and Mobley (2003), which provides a polynomial parameterization of the deep water reflectance R rs∞ (λ) (units: sr −1 ) in terms of the dimensionless backscattering albedo x (λ): with coefficients p 1 = 0.0512 sr −1 , p 2 = 4.6659 sr −1 , p 3 = −7.8387 sr −1 , p 4 = 5.4571 sr −1 , p 5 = 0.1098 sr −1 , p 6 = −0.0044 s m −1 , and p 7 = 0.4021 sr −1 . θ s is the zenith angle of the sun under the water surface and θ v is the sensor viewing angle also below the water surface, and u is the mean wind speed. We used a mean value over the world's oceans of u = 6.9 m s −1 , θ s = 29 • , typical for summer at temperate latitudes (40-60 • N). The model is valid even in turbid water with large backscattering albedo, which is given by: with total absorption coefficient, a (λ)(in m −1 ), composed as follows a (λ) = a w (λ) + a cc (λ) + a cdom (λ) + a ph (λ) where a w (λ) , a cc (λ) , a cdom (λ), and a ph (λ) are the respective contributions from pure seawater, coccolithophore cores, colored dissolved organic matter (CDOM), and phytoplankton other than coccolithophores. To focus on the effect of coccolithophores, we assumed a w (λ) , a cdom (λ), and a ph (λ) to be constant. We fixed the chlorophyll-a concentration from other phytoplankton to a value of 0.6 mg m −3 and set a cdom = 0.02 at a wavelength of 0.443 µm with a spectral slope of 0.0172 nm −1 , corresponding to average values for the Atlantic (Babin et al., 2003). See Supplementary Materials for more details on the calculation of a(λ). The total backscattering coefficient b b (λ) is partitioned as follows: where b bw (λ)is the backscattering coefficient of pure water and b bpic (λ) is the backscattering coefficient of the coccolithcoccosphere mixture as in Equation (1). We note that we have assumed negligible contribution to backscattering from the organic part of the coccolithophore as the refractive index of the organic part of the coccolithophore is estimated to be 1.06 relative to water (Aas, 1996), much smaller than the refractive index of its calcite coccosphere, which is 1.20 relative to water. Further details of the model can be found in Supplementary Information section 1.3. As derived in Supplementary Information section 1.3, the coccolith core absorption term is given by where N cc is the number concentration of coccospheres (in m −3 ) and a * cc (λ) is the coccosphere-specific absorption coefficient (in m −2 per coccosphere) resulting from an integral over the size distribution of the cores: We note that if we neglect the backscatter loss from the front surfaces which we can do to first order, the absorption by the cores is independent of whether they are coated or not. We can thus write: where F nc is the proportion of naked cores and F sc is the fraction of naked cores that remain in the surface after their liths have been shed (F sc = number of surviving naked cores /N cc ). The first term represents the absorption from the naked cores left in the surface layer after lith shedding, while the second term represents the absorption from the coated cores in the bloom. The backscattering term results from Equation (54): Even though cell densities in bloom conditions for E. huxleyi are somewhat arbitrary, a minimal value of 10 9 cells m −3 has been suggested (Tyrrell and Merico, 2004), with cell concentration exceeding 10 11 cells m −3 in the most intense blooms observed in Norwegian fjords (Tyrrell and Merico, 2004). Therefore, a realistic range for initial coccolithophore concentrations would be between N cc = 10 9 and 2 × 10 10 cells m −3 in open ocean waters, which we will use to simulate changes in remote sensing reflectance.

Spectral Backscattering Efficiency for a Single Coccolith as a Function of Size, Shape, and Morphology
The refractive index of calcite does not vary significantly over the wavelength range and can be represented by an average value of 1.2, which gives ω brg and ω bsm (Fournier and Neukermans, 2017), so that Q bbl is only a function of λ w /a o , as shown in Figure 4. For a coccolith with semi-major axis length a o , in the long wavelength limit the entire surface of the particle appears perfectly smooth since asperities or holes will be smaller than λ w /4, and Q bbl = ω bsm 1 − RF rg . As the wavelength gets shorter, parts of the surface will start to appear as rough and Q bbl will increase rapidly to reach a maximum value of Q bbl = ω brg F rg + ω bsm 1 − F rg when all asperities and holes are larger than λ/4. This is very similar to the result obtained when modeling circular coccoliths with the morphology described by Zhai et al. (2013). The backscattering efficiency for an elliptical coccolith with various aspect ratios r 0 = a o /b o is shown in Figure 4A. The change seen in the onset wavelength is due to the reduction in the maximum gap width as the overall perimeter of the lith is reduced because the axis ratio increases. The accompanying change in the maximum value of the backscattering efficiency is due to the increase in the maximum fraction F rg of rough area to total area because of the increasingly elliptical shape. The effect of the gap shape factor r sh (defined in Equation 9) is shown in Figure 4B for an elliptical axis ratio r 0 of 1.2. This factor is a significant determinant of the spectral shape of light backscattering. In the limiting case of a constant gap width (r sh ≈ 1) the complete change from smooth to rough backscatter will occur in a very narrow wavelength range. As the wedge gap narrows to a triangular shape (r sh ≈ 0) the change will become more gradual. As this triangular shape is the one observed by Young et al. (2014), we will use it in the rest of this paper. Figure 5 shows the spectral shape of backscattering efficiency for various coccolith size distributions modeled as a second order shifted Gamma function with five different mean semi-major axes µ and five standard deviations σ as defined in Equation (27). Not surprisingly, the sharp spectral features observed for a single coccolith particle in Figure 4 are smoothened when averaged over the size distribution. Figure 5 shows that the spectral shape of backscattering is largely controlled by variations in mean size of the coccolith population, whereas the standard deviation is a secondary factor merely broadening the transition zone from smooth to rough scattering as the standard deviation increases.

Spectral Backscattering Efficiency for Populations of Elliptical Coccoliths
Differences in σ bbl (λ) between the first order gamma function used in our original work the second-order gamma models for coccolith size distributions are very small because the Gamma functions have the same mean and standard  deviation. However, in what follows, we will use the second-order Gamma function to represent the coccolith size distribution as it matches the experimental size distributions more closely (see Supplementary Information section 1.2).
The changes in the spectral shape of light backscattering efficiency as the mean lith size changes are due to the fact that the transition from smooth to rough backscattering occurs at longer wavelengths as the size of the gap increases, which corresponds to an increase in coccolith size for a fixed coccolith morphology. The effect is shown in Figure 5A for a change in mean coccolith major axis of 2.0 to 3.6 µm, which corresponds to the range observed in the lab and the field.

Differences in Spectral Light
Backscattering Between Free and Attached Coccoliths Figure 6A depicts variations in the spectral behavior of the ratio σ bbNl /σ bbs for a coccolithophore with a single layer of liths and core radius r c ranging between 2.75 and 3.2 µm which matches the ranges observed by Young et al. (2014). At wavelengths shorter than 0.5 µm, backscattering from free liths is 40-60% higher than backscattering from coated spheres and about 30% higher at larger wavelengths. Increasing the number of coccolith layers reduces the region where the core is absorbing and leads to smaller differences in backscattering cross-section between free and attached liths as can be seen in Figure 6B. We used a core radius of 3.0 µm which corresponds to a mean lith major axis value of 3.3 µm (Zhai et al., 2013), close to those measured by Young et al. (2014). For four layers of coccoliths, the total number of coccoliths on the coccosphere is 40. For that many layers the coated coccolithophores are nearly spectrally indistinguishable from the free liths and backscattering from free liths is at most 35% higher than backscattering from coated spheres.
Note that for a single layer of interlocked liths the first surface could simply be completely covered. If there are multiple layers however, all the layers must have about the same number of coccoliths per layer. The interlocked liths model leads to substantially smaller number of liths for a given number of layers and reduces further the difference between the backscattering from the sum of the free coccoliths and from the fully coated coccolithophore cores.

Changes in Spectral Light Backscattering
During E. huxleyi Coccolith Shedding Event Figure 7 shows changes in the magnitude and spectral shape of the backscattering cross section σ bbtot (λ) for a coccolithophore with a core radius of 3.7 µm, shedding four layers of liths with mean a o of 1.7 µm and standard deviation of 0.15 µm, obtained through integration of expressions in Equations (55) and (57). Note that the particulate backscattering coefficient b bp (λ) can be obtained by multiplication of σ bbtot (λ) with the initial number concentration of coated cores, N cc . The size of the coccosphere (which also fixes the size of the liths) and number of coccolith layers chosen for this simulation corresponds to a coating of the coccosphere with 36 liths ( Table 1). This gives a coccolith-tocell ratio of 36, on the lower range of what has been observed for coccolith-to-cell ratios in natural and laboratory E. huxleyi blooms, which vary between 45 and 400 liths per cell (e.g., Balch et al., 1991Balch et al., , 1993Gordon et al., 2009). Potential causes for a lower coccolith-to-cell ratio than observed in laboratory experiments and in the field are: (i) E. huxleyi's continued calcification during the coccolith shedding phase, (ii) the presence of detached coccoliths before the shedding begins , and (iii) the preferential sinking of coccospheres out of the surface layer due to their larger size, but see Monteiro et al. (2016).
The spectral shape of b bp (λ) is commonly described by the relationship (Gordon and Morel, 1983): where λ o is a reference wavelength and ξ represents the dimensionless spectral slope parameter of b bp (λ). Calculated over the 0.4-0.7 µm wavelength range, ξ changes from −0.77 to −0.94 for coated spheres to free liths ( Figure 7B). Even though not strictly quantitatively comparable, our results are consistent FIGURE 6 | Ratio of free coccolith to coated coccolithophore backscattering cross-section σ bbNl /σ bbs , for (A) single-layer coated coccospheres of varying core radii and (B) coccospheres with varying number of coccolith layers and a core radius of 3.0 µm.
Frontiers in Marine Science | www.frontiersin.org FIGURE 7 | Changes in the total spectral backscattering cross-section σ bbtot magnitude (A) and spectral shape (B) as a coccolithophore with a core radius r c of 3.7 µm, a lith semi-major axis a 0 of 1.7 µm, and 4 layers of liths sheds its liths.
with multispectral b bp (λ) measurements in E. huxelyi cultures who reported a steepening of the spectral slope of b bp (λ) for suspensions of free liths compared to coccospheres (Voss et al., 1998), from ξ = −1.2 to ξ = −1.4. We note that the experimental error on the slopes calculated by Voss et al. (1998) were large and that the morphology of E. huxelyi was not determined. A change to the shape or size of the openings will lead to a change of slope as shown in Figure 4 for a single lith.

Ocean Color Remote Sensing Reflectance for E. huxleyi Populations of Different Size Distribution and Bloom Stage
Implementing our optical model for light backscattering from E. huxleyi coccoliths and coccospheres into the analytical radiance model presented in section Algebraic Radiance Model, we can investigate the variability in ocean color remote sensing reflectance for blooms of different magnitudes, with different lith size distributions, and during different stages. The variability in magnitude and spectral shape of the abovewater remote sensing reflectance R rs (λ) (in sr −1 ) during a lith shedding event of E. huxleyi blooms with initial coccolithophore concentrations ranging between N cc = 10 9 and 2 × 10 10 cells m −3 is depicted in Figure 8. We assumed a coccolith size distribution with a mean a 0 of 1.7 µm and standard deviation of 0.15 µm and coccospheres with four layers of liths, totaling 36 coccoliths per coccosphere for all our simulations, unless noted otherwise. Figures 8A,B shows the magnitude and spectral shape of R rs (λ) at the initial stage with all liths attached (F nc = 0) and Figures 8C,D shows the final stage with all liths shed (F nc = 1) and no living cores remaining (F sc = 0). Before the onset of lith shedding, the absorption features of the cores clearly influence R rs (λ) by reducing reflectance in blue wavelengths (∼ 0.4 − 0.47 µm), resulting in peak R rs (λ) at green wavelengths (∼ 0.50 − 0.55 µm). When the liths are freed and no cores are left to absorb light, R rs (λ) increases by a factor of up to four and the peak in R rs (λ) shifts toward blue-turquoise wavelengths (∼ 0.40 − 0.49 µm), where R rs (λ) shows little spectral variation. This gives a bright milky turquoise coloring of the ocean, typically associated with E. huxleyi blooms. Our modeled R rs (λ) spectra are also quantitatively consistent with multispectral ocean color satellite observations of R rs (λ) in coccolithophore blooms which have been shown to peak at 0.49 µm (Moore et al., 2012).
Next we investigate the influence of changes in the fraction of naked cores (F nc ) or surviving cores (F sc ) on R rs (λ) during different stages of a bloom with initial E. huxleyi coccolithophore concentration of 10 9 cells m −3 . Figures 9A,B shows how R rs (λ) changes when the coated cores shed their liths and they die after shedding, i.e., F sc = 0. Even though there are significant changes in the magnitude of R rs (λ), the spectral variability is clearly small, and appears to be too small to be detectable from satellite ocean color observations which typically have uncertainties around 5% at visible wavelengths. The change in spectral shape of R rs (λ) due to lith shedding will, however, be more pronounced for larger cell concentrations. Figures 9C,D illustrates variability in R rs (λ) due to changes in the fraction of cores that remain in the surface after lith shedding (F sc from 0 to 1 when F nc = 1). Such changes result in much stronger spectral variability in R rs (λ) by virtue of their absorption features (Figures 9C,D).
Lastly, Figure 10 illustrates how changes in the mean size of free liths impact the magnitude and spectral shape of R rs (λ). The magnitude of R rs (λ) increases with increasing mean size mostly due to increases in the geometric cross section. The spectral shape of R rs (λ) varies significantly over the wavelength range considered: smaller-sized coccoliths produce stronger relative reflection at shorter wavelengths, which offers the potential for ocean color remote sensing of coccolith size distribution.

DISCUSSION
Our model results indicate that light backscattering from detached E. huxleyi coccoliths is at most twice as large as the backscattering from the same liths when arranged in a coccosphere. This result is consistent with the suggestion of Gordon et al. (2009), but may come as a surprise as it has been commonly assumed that the backscattering from coccoliths would be orders of magnitude higher than for coccospheres. The implication of this result is that the calcite-specific backscattering coefficient, b bp /PIC, for E. huxleyi coccoliths of a given size and morphology varies by at most a factor of two between coated cores and free liths. For E. huxleyi coccolith populations of a given size and morphology, light backscattering is thus a good proxy for calcite mass concentration, regardless of whether coccoliths are attached to or freed from the core. The difference between scattering from free and attached liths is smallest in regions of the spectrum where absorption of the core is minimal, i.e., between 0.54 and 0.66 µm and wavelengths larger than 0.7 µm. Another implication is that the order of magnitude increase in backscattering frequently observed at the final bloom stages is likely mostly due to the excess production of coccoliths and only a small part of this increase is due to the shedding of the attached liths.
We have demonstrated that the magnitude and spectral shape of light backscattering by coccoliths are strongly sensitive to coccolith size and morphology, which suggests that information on the size and morphology are potentially amenable from hyperspectral backscattering measurements. In what follows, we describe a potential avenue for estimating coccolith size from b bp (λ). The spectral shape of b bp (λ) from a single coccolith is characterized by the onset of a sharp increase at a definite wavelength due to a transition from smooth to rough scattering. This onset wavelength in air, λ onset , is determined by the gap width at the outer perimeter from Equation (10) as follows: where N is the number of gaps, r pg is the pillar to gap ratio defined in Equation (9), n w is the refractive index of seawater, and P n is the non-dimensional form of the perimeter of the lith defined in Equation (S13). For a given lith morphology, which fixes the values of N, r pg , a r /a o , and b r /b o , the lith size a 0 can thus be estimated. Equation (67) expresses the condition that at the transition from smooth to rough backscattering, the gap width is equal to a quarter of the wavelength in water. As the wavelength gets shorter the backscattering efficiency increases from a value of ω bsm to a maximum value of: where F rg is the fraction of the surface of the lith which is composed of gaps and pillars and which obviously depends on the assumed morphology (see Equation 7), and ω brg and ω bsm are respectively the rough and smooth reflectivities of the lith FIGURE 9 | Spectral above-water remote sensing reflectance R rs (λ) in absolute values of sr −1 (A,C) and normalized at 0.55 µm (B,D) for various fractions of naked cores (F nc = 0 to 1, top row), and various surviving core fractions (F sc = 0 to 1, bottom row). We used a size distribution of free liths of mean a 0 1.7 µm and standard deviation of 0.15 µm, an initial coccolithophore cell concentration of 10 9 m −3 and four layers of liths. surfaces and only depend on the index of refraction of the calcite material of the lith, which is 1.2.
In the simple case of a population of coccoliths of fixed morphology, the parameters that directly control the spectral shape of b bp (λ) are and F rg . Our results suggest that the mean lith size of the population may be obtained from the spectral shape of b bp (λ) assuming a functional form for the size distribution of liths expressed in terms of a mean particle size µ and standard deviation σ . To help solve for the size distribution parameters, it would be advantageous if a relationship between µ and σ could be established. From proportionality arguments one would expect a linear relationship to apply to first order. Such a fit applied to the data of Young et al. (2014) gives σ ≈ 0.11µ . However, because of the small range of sizes and the significant spread of σ the R 2 coefficient of determination is only 0.25. More particle size distribution measurements over larger size ranges will be needed to improve solving for mean population size.
All the above relies on an assumed lith morphology, which is known to vary among morphotypes of E. huxleyi (Young et al., 2003;Cook et al., 2011;Hagino et al., 2011). The coccolith morphology and size range assumed in this paper are based on the work of Young et al. (2014) for E. huxleyi of morphotype A. Optical models for different morphotypes of E. huxleyi could be set up for example by adjusting the proportionalities in Equations (6) and (9). The corresponding inversion algorithms to estimate lith size from b bp (λ) would thus also be morphotypespecific. Of primary importance in these efforts is the relationship between the gap width and the lith size, as this forms the basis of the size estimation. As pointed out earlier, an assessment of the relationships between statistical parameters of the size distribution, such as the standard deviation and mean, would also be desirable.
Our model further illustrates that when a population of E. huxleyi sheds its liths, the spectral shape and magnitude of b bp (λ) change surprisingly little. This is due to the quasiequality of the backscattering cross-sections of the same amount of liths when free or arranged in a coccosphere. The subtle changes in b bp (λ) are due to the increase in backscattering by the free liths that are no longer obscured by the core absorption. Nevertheless, our results suggest that estimation of the fraction of cores that have shed their liths (F nc ) appears to be feasible if accurate measurements of the spectral slope of b bp (λ) can be obtained. This information may be exploited to develop approaches to determine the proportional contribution of coccoliths and coccolithophores to total calcite, which may be of interest in biogeochemical and ecological studies of coccolithophores. While coccolithophores are marine microalgae which contribute to primary production, calcification (Monteiro et al., 2016), and the production of dimethylsulfoniopropionate (Matrai and Keller, 1993), coccoliths are dead material thought to enhance the flux of organic carbon to the deep ocean by providing ballast material (Armstrong et al., 2001;Riebesell et al., 2016).
We have also modeled the remote sensing reflectance R rs (λ) during various stages of hypothetical E. huxleyi blooms in which coccoliths are shed and the naked cores either die or survive. This allows us to investigate which characteristics of E. huxleyi bloom state are potentially resolvable from hyperspectral ocean color remote sensing. As expected from the subtlety of changes in spectral b bp (λ) when coccolithophores shed their liths and naked cores die, changes in spectral R rs (λ) are similarly small and therefore unlikely to be detectable from space. In a population of free liths at the end of a shedding event, the fraction of surviving cores had a much more pronounced effect on R rs (λ) by virtue of their absorption. However, similar effects on R rs (λ) can be expected from changes in the absorption of phytoplankton other than coccolithophores, which is assumed constant in our simulations. Therefore, the potential to estimate naked E. huxleyi cells depends on how different their absorption features are from other types of phytoplankton.
Our modeled R rs (λ) were produced assuming a constant background of absorption by CDOM and chlorophyll-a from phytoplankton other than E. huxleyi, as well as negligible contribution to backscattering from particles other than E. huxleyi coccoliths and coccospheres. This is obviously an overly simplified representation of real bloom conditions. The assumption of constant non-zero CDOM absorption was essential to obtain the spectral shape of R rs (λ) typically associated with E. huxleyi blooms, i.e., one that produces a turquoise hue with peak reflectance at 0.49 µm (Moore et al., 2012). This is because CDOM absorbs strongly in violet-blue wavelengths. If CDOM absorption is not included in the background, R rs (λ) would peak at violet wavelengths and the magnitude of R rs (λ) would be much more intense than observed from ocean color satellites and in the field. The impact of CDOM absorption on reflectance from a suspension of 3.6 × 10 10 m −3 free coccoliths and no absorption due to chlorophyll-a is demonstrated in Figure 11A, whereas Figure 11B also includes the absorption effect of 10 9 m −3 E. huxleyi cores.
In coccolithophore blooms characterized by high concentrations of PIC (>0.003 mol C m −3 ) the ocean color FIGURE 11 | Spectral remote sensing reflectance R rs (λ) (in sr −1 ) from the free liths shed by N cc = 10 9 E. huxleyi cells m −3 in waters with a ph (λ) =0 and varying CDOM absorption, when (A) all cores have died, F sc = 0 and when (B) all cores survived, F sc = 1. The coccolith size distribution had a mean of 1.7 µm and standard deviation of 0.15 µm and the coccospheres were coated with four layers of liths. remote sensing algorithm of  is used to derive PIC. Because bands at shorter wavelengths often saturate in these bloom conditions, the algorithm uses R rs (λ) in three bands in the red and near-infrared (0.670, 0.765, and 0.865 µm) to estimate b b from R rs at 0.670 µm assuming R rs = 0 in the 0.765 and 0.865 µm bands to perform the aerosol scattering correction. From the range of coccolith concentrations used in Figure 8C we find R rs at 0.765 and 0.865 µm much lower than R rs at 0.670 µm, thereby providing support for the atmospheric correction approach of . We suggest that our model may be used to improve on the atmospheric correction approach of .
We have assumed throughout this work that the detached coccoliths are randomly oriented with respect to incident light, which allows us to apply the theorem of Cauchy (1850) to calculate the average geometric cross-section via Equation 3. However, it has recently been shown that large phytoplankton cells and marine aggregates have a tendency to orient horizontally in natural waters (Nayak et al., 2018). If coccoliths, even though much smaller in size, have just a slight tendency to orient horizontally with respect to incident light, the influence on backscattering will be significant, because light is scattered off a larger surface. As a consequence the difference in backscattering between oriented coccoliths and coccolithophores will be much larger. We can estimate the limits of this effect by taking the ratio of the randomly oriented liths surface to the surface at normal incidence.
R ori ≈ πab 2π ab+2π (a+b)τl 4 = 2 1 + (1 + r) τ l a (69) For completely oriented very thin liths the cross-section is therefore approximately twice that of randomly oriented liths. To distinguish this orientation effect in natural waters, new observational approaches such as circular polarization will be needed to study preferential particle orientation.
Besides the assumption of random orientation of the coccoliths and coccolithophores, our approach is also based on the assumption that the unity-normalized angular distribution of reflection from any randomly oriented convex body is independent of shape and size (van de Hulst, 1957), implicit in Equation (3). When backscattering is dominated by reflection, such as for high refractive index coccoliths and coccolithophores, this implies that the scattering efficiency Q bb is independent of shape and size and only dependent on the relative amount of rough and smooth surfaces comprising the particle. This angular constancy is the physical underpinning of our analytical light backscattering model and also allows us to carry out the averaging of Q bb over a size distribution. The fact that our modeled backscattering cross-sections for circular disk-like coccoliths compared successfully with exact scattering solutions from the discrete dipole approximation (Fournier and Neukermans, 2017) provides general support to our assumptions. We expect the validity of the angular constancy assumption to break down as the particle shapes become more complex and self-shadowing (non-convexity) becomes significant, which can be investigated using exact light scattering codes such as discrete dipole approximation, T-Matrix, or Finite Difference Time Domain (FDTD) approaches (e.g., Gordon and Du, 2001;Hedley, 2012;Zhai et al., 2013).
A thorough validation of our light scattering model for E. huxleyi requires simultaneous measurements of E. huxleyi coccolith morphology, size distribution, and spectral light backscattering, which can be obtained in laboratory cultures. Further validation work will also entail conjunct in situ observations of (hyperspectral) remote sensing reflectance, inherent optical properties, E. huxleyi coccoliths and coccolithophore size distribution and morphology, and characterization of other optically active constituents during E. huxleyi blooms. Such measurements in the laboratory and in natural waters are essential to the evaluation of the current model and will inform on potential modeling improvements.