METHODS article

Front. Mar. Sci., 07 May 2018
Sec. Ocean Observation

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

  • 1Laboratoire d'Océanographie de Villefranche, Sorbonne Université, Centre National de la Recherche Scientifique, Paris, France
  • 2DRDC, Valcartier Research Laboratory, Valcartier, QC, Canada

In this study we develop an analytical model for spectral backscattering and ocean color 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 coloring 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.


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 (Gordon et al., 2001; 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 (Gordon et al., 2001; 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; Neukermans et al., 2018).


Figure 1. Scanning Electron Micrograph of Emiliania huxleyi (Dr. Jeremy Young, University College London, London, with permission).


Table 1. List of abbreviations, symbols, definitions, and units.

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, 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.


The light backscattering coefficient of calcite particles suspended in seawater, bbpic (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 bbl and bbc, respectively:

bbpic=bbl+bbc    (1)

For populations of identical coccoliths and coccospheres we can write:

bbpic=Nlσbbl+Nccσbbs    (2)

where Nl and Ncc are the number concentrations of coccoliths and coccospheres (in m−3), respectively, and σbbl and σbbs are the corresponding backscattering cross-sections for individual coccoliths and coccospheres (in m2). The backscattering cross-section is defined as the product of the geometric cross-section of the calcite particle σg (in m2) and the dimensionless efficiency factor for backscattering of the particle Qbb so that:

σbbl=Qbblσglandσbbs=Qbbsσgs    (3)

where the efficiency factor Qbb 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 Qbb 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 Qbb.

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:

Qbbl=ωbsm(1RFrg)+ωbrgRFrg    (4)

where Frg is the maximum fraction of the projected area of the particle that can become rough, RFrg 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 Qbbl 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 wr 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 ww 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 wt that represents the boundaries of the inner tube joining the distal and proximal sheets.


Figure 2. Schematic representation of the morphology of an E. huxleyi coccolith as measured by Young et al. (2014). The inset describes details of the wedge structure of the openings.

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:

σgl=(14)[2π(ab)+τl Pe(a,b)]    (5)

where τl is the thickness of the elliptical disk and Pe(a, b) its perimeter with equations given in Supplementary Information section 1.1.

We can represent the proportional relations between the various ellipse parameters (ao, ar, aw, at) and (bo, br, bw, bt) as follows (Figure 2):

arao=[1wra0] and brbo=[1(aobo)wrao]awao=[1(wr+ww)a0] and bwbo=[1(aobo)(wr+ww)ao]atao=[1(wr+ww+wt)a0] and btbo=[1(aobo)(wr+ww+wt)ao]    (6)

For the coccolith in Figure 1, with major axis in μm of ao = 1.8688, we have the following proportionalities:


The maximum fraction of the surface area of the coccolith that can become rough Frg lies between ellipses with parameters (ar, aw) and (br, bw) so we can write:

Frg=[araobrboawaobwbo]    (7)

which gives Frg = 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 ar and aw can be represented as a set of gaps with maximum widths Δsgmax and minimum pillar width Δspmin respectively along the elliptical boundary at ar and minimum widths Δsgmin and maximum pillar width Δspmax. The gap width narrows down to Δsgmin on the elliptical boundary at aw while the pillar width grows to fill the empty space. Given Np pillars around the circumference of the ellipse at ar we have:

Np(Δsgmax+Δspmin)=Pe(ar,br)Np(Δsgmin+Δspmax)=Pe(aw,bw)    (8)

Defining the ratio of the minimum pillar width to the maximum gap width, rpg, and the ratio of minimum to maximum gap width, rsh, as:

rpg=ΔspminΔsgmax    and  rsh=ΔsgminΔsgmax    (9)

we obtain the following expression for Δsgmax

Δsgmax=Pe(ar,br)Np(1+rpg)=2πa0 Pn(arao,brao)Np(1+rpg)    (10)

where Pn is the dimensionless ellipse perimeter given in Equation (S13). As we proceed from the ellipse at ar to the ellipse at aw the gap Δsg narrows from Δsgmax to Δsgmin as w goes from 0 to ww

Δsg=Δsgmax[1(1rsh)www]    (11)

To determine where the backscattering becomes rough we need to compute at what value of a the gap Δsg becomes equal to 1/4 of λw (Gordon, 2006; Fournier and Neukermans, 2017). We define wc as the critical width where the gap becomes equal to λw/4, then Equation (11) gives:

λw4=Δsgmax[1(1rsh)wcww]    (12)

The corresponding critical semi-major axis ac and semi-minor axis bc are given by:

ac=ao(wr+wc) or acao=[1(wr+wc)ao]    (13)
bc=bo(wr+wc) or bcbo=[1(aobo)(wr+wc)ao]    (14)

As long as the gap Δsg is smaller than Δsgmax and larger than λw/4 we can compute the critical width as:

wc=ww(1λw/4Δsgmax)(1rsh)   for   λw(4Δsg)1    (15)

When Δsg is smaller than λw4 we must set wc to 0 and on the high side as the wavelength gets smaller we must limit wc to a maximum value of ww:

wc=0 for λw/4Δsg>1 and wc=ww for (1λw/4Δsgmax)(1rsh)>1    (16)

We can then proceed to compute the ratio wcao which is required to evaluate the fraction of rough surface of the coccolith using Equations (15) and (10).

wcao=wwao(11rsh)(1λw/4Δsgmax)      =wwao(11rsh)[1(λw4πao)Np(1+rpg)Pn(arao,brao)]    (17)

Substituting this expression in Equations (13) and (14) we obtain directly the fraction of the elliptical surface of the coccolith that is rough:

RFrg(ao,λ)=πarbrπacbcπaobo=[araobrboacaobcbo]    (18)

The only variable determining the behavior of RFrg is the ratio of the wavelength in water λw to the major axis ao 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:

Pe(aw,bw)Pe(ar,br)=11+rpg(rsh+ΔspmaxΔsgmax)=11+rpg(rsh+fpbrpg)    (19)


fpb=Δspmax/Δspmin    (20)

is the ratio of maximum to minimum pillar width. For triangularly shaped wedges such as those measured by Young et al. (2014) Δsgmin = 0 = rsh and using Equation (S13) for a non-dimensional version of the perimeter ratio, Pn, we find:

Pe(aw,bw)Pe(ar,br)=Pn(aw/ao,bw/a0)Pn(ar/ao,br/a0)    (21)

Equation (19) then simplifies to:

Pn(aw/ao,bw/a0)Pn(ar/ao,br/a0)=fpb rpg1+rpg    (22)

For the coccolith morphology studied by Young et al. (2014) and shown in Figure 1 we have rpg = 0.666 and fpb = 1.602 which means that the base of the pillar Δspmax is almost as wide as the wedge opening at the top Δsgmax. 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):

σgl=(14)[2π(aobo)+τl πao Pn(aoao,boao)]       =(πaobo2)[1+τl2bo Pn(1,boao)]    (23)

with Pn 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:

τn=12(aob0)τlao    (24)

which gives

σgl=(πaobo2)[1+τn Pn(1,boao)]


σgl=(πao22)(boao)[1+τn Pn(1,boao)]    (25)

In our previous work (Fournier and Neukermans, 2017) we assumed a first-order shifted gamma size distribution P(ao) 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:

P(ao)=4β3/2π(aoaom)2eβ(aoaom)2aoaom    (26)

where aom 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:

β=(3π8)2π σ2, aom=μ  σ(3π8)8, apeak=μ  σ(3π8)2 π    (27)

where apeak 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:

σgl=π2(aobo)[1+τn Pn(1,boao)]aomao2 P(ao) dao    (28)


σgl=π2(aobo)[1+τn Pn(1,boao)][aom2+4aom1π β+32 β]    (29)

The mean backscattering cross-section and backscattering efficiency factor for a population of elliptical coccoliths are then given by, respectively:

σbbl(λ)=aomQbbl(λao)σgl(ao) P(ao) daoQbbl(λ)=aomQbbl(λao) P(ao) dao    (30)

Note that P(ao) can now be expressed as a pure function of measured mean μ and standard deviation σ by using the relations in Equation (27):

P(μ,σ,ao)=4β(σ)32π(aoaom(μ,σ))2eβ(aoaom(μ,σ))2aoaom    (31)

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 randomly oriented convex shape is equal to ¼ of its surface area (Cauchy, 1850) we find for a coccosphere of radius rs that

σgs=(14)4πrs2=πrs2    (32)

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).

For a randomly oriented elliptical disk of semi-major axis a, semi-minor axis b, aspect ratio r = a/b, and thickness τd, the mean geometric cross-section is:

σgl=(14)[2π(ab)+π(a+b)τd]=πa22r[1+(1+r)2rτda]       =πa22r[1+(1+r)2rfd]    (33)

where fd = τ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 Pe(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:

Nπab=Nπa2r=4πrc2    (34)

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., Qbbs = 2Qbbl, and the total backscattering cross-section of the coccosphere is given by:

σbbs=2Qbblπrs2=Qbbl[Nπa22r ]    (35)

Using Equations (33) and (34), the backscattering cross-section for N disks, σbbNl, is given by:

σbbNl=Qbbl[Nπa2(1+(1+r)2rfd)2r]             =Qbbl[4πrs2(1+(1+r)2rfd)2]    (36)

which can be simplified to:

σbbNl=Qbbl[2πrs2(1+(1+r)2rfd) ]    (37)

From Equations (35) and (37), we thus find that for a non-absorbing 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:

σbbNlσbbs=1+(1+r)2rfd    (38)

This expression is near unity for thin disks (fd ≈ 0).

To account for the multi-layered structure of the coccosphere, let Ovl be the number of layers. Adding (Ovl − 1) additional layers of coccoliths attached by their proximal sheets to the distal sheets of the coccoliths on the layer beneath them we have:

rs2=(rc+(Ovl1) Δrs)2    (39)

where Δrs 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 Ovl number of layers:

σbbNlσbbs=2πrs2(1+(1+r)2rfd)2πrs2=1+(1+r)2rfd    (40)

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):

σbbNlσbbs=2(1+(1+r)2rfd)[1+(1Qas(k,2rc)rc2rs2)]    (41)

Where Qas(k, 2rc) 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):

Qas(k,2rc)=2[12+eδ2 rcδ2 rc+(eδ2 rc1)(δ2 rc)2] with δ=4k 2πλ    (42)

To evaluate Qas(k, 2rc) we use the spectral coccolithophore core absorption k values obtained by Stramski et al. (2001) and Neeley et al. (2015). Using Equation (39), we can now rewrite Equation (41) as a function of Ovl:

σbbNlσbbs=2(1+(1+r)2rfd)[1+(1Qas(k,2rc)1(1+OvlΔrsrc)2)]    (43)
σbbNlσbbs=2(1+(1+r)2rfd)[2Qas(k,2rc)1(1+OvlΔrsrc)2]    (44)

So in the case of thin liths (fd ≈ 0) and layers (rs ≈ 0) and a fully absorbing core (Qas(k, 2rc) ≈ 1) this would imply that the scattering from the free liths would be twice as much as the scattering from the coated coccosphere:

σbbNlσbbs2    (45)

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 chlorophyll-a 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 aw, we find the following equality:

Nil π a02r(1awbwa0b0)=Nil π a02r(1ft2)=4πrc2    (46)

where a0 is the semi-major axis of the outer part of the coccolith, rc is the radius of the coccolithophore core, and Nil 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, Nsl, we find the ratio of the number of liths per layer for the two models:

NilNsl=(1ft2)    (47)

Expressing the exclusion factor ft2 in terms of the geometry of the coccolith we have:

ft2=awbwa0b0=[1(wr+ww)a0][1r(wr+ww)a0]    (48)

For liths with the morphology measured by Young et al. (2014) we have


and therefore

NilNsl=0.6365    (49)

Measurements on cultured E. huxleyi cells give typical values for the mean core radius rc of 2.2 μm and the mean semi-major axis a0 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 a0/rc 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 aeq:


we find that aeq/rc 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/rc =0.07), the separation between the sheets (dh/rc =0.18), and the tube width (at/rc =0.18).

Solving Equation (34) for Nsl we find:

Nsl=4 r rs2a02=4× 1.2×(2.21.25)2=14.8715

And using Equation (49) we have


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, Δrs, should be 1.53 μm. For l layers this gives

Nsl=4 r (rc+l Δrs)2a02    (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.


Table 2. Comparison of the number of coccoliths, Nl, attached to a coccosphere of radius 3.7 μm for 1–4 layers for stacked and interlocked liths of semi-major axis length 1.7 μm and as observed on 3D imagery of coccospheres by Hoffmann et al. (2015).

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:

2πrc2(1ft2)Ovl    (51)

We thus find the following ratio for a non-absorbing core:

σbbNlσbbs=Qbbl 2πrc2(1ft2)Ovl (1+(1+r)2rfd)Qbbl 2πrc2(1ft2)Ovl =1+(1+r)2rfd    (52)

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:

σbbNlσbbs=2(1+(1+r)2rfd)[2Qas(k,2rc)1(1+OvlΔrsrc)2]    (53)

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):

fd=14 , Δrsrc=15 , r=1.2

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 (Balch et al., 1993). 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 (Balch et al., 1993) 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 Nltot liths attached. This suspension then gradually transforms into a suspension comprised of free liths with a certain fraction of surviving absorbing cores left, Fsc, which varies between 0 and 1.

If we define Fnc as the ratio of cores that have shed their liths to the total initial number of coated cores (Fnc = number of naked cores/Ncc), 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:

σbbtot=FncσbbNl+(1Fnc)σbbs    (54)


σbbNl=Qbbl2πrc2(1ft2)Ovl(1+(1+r)2rfd)    (55)

the backscattering cross-section of the free liths when detached from a single coccolithophore, and

σbbs=σbbNl[2Qas(k,2rc)1(1+OvlΔrsrc)2]2(1+(1+r)2rfd)    (56)
σbbs=Qbblπrc2(1ft2)Ovl[2Qas(k,2rc)1(1+OvlΔrsrc)2]    (57)

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 Rrs(λ) (units: sr−1) in terms of the dimensionless backscattering albedo x(λ):

Rrs(λ)=p1x(λ)(1+p2x(λ)+p3x(λ)2                          + p4x(λ)3)(1+p51cosθs)(1+p71cosθv)(1+p6u)    (58)

with coefficients p1 = 0.0512 sr−1, p2 = 4.6659 sr−1, p3 = −7.8387 sr−1, p4 = 5.4571 sr−1, p5 = 0.1098 sr−1, p6 = −0.0044 s m−1, and p7 = 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:

x(λ)=bb(λ)bb(λ)+a(λ)    (59)

with total absorption coefficient, a(λ)(in m−1), composed as follows

a(λ)=aw(λ)+acc(λ)+acdom(λ)+aph(λ)    (60)

where aw(λ), acc(λ), acdom(λ), and aph(λ) 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 aw(λ), acdom(λ), and aph(λ) to be constant. We fixed the chlorophyll-a concentration from other phytoplankton to a value of 0.6 mg m−3 and set acdom = 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 bb(λ) is partitioned as follows:

bb(λ)=bbw(λ)+bbpic(λ)    (61)

where bbw(λ)is the backscattering coefficient of pure water and bbpic(λ) is the backscattering coefficient of the coccolith-coccosphere 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

acc(λ)=Nccacc*(λ)    (62)

where Ncc is the number concentration of coccospheres (in m−3) and acc*(λ) is the coccosphere-specific absorption coefficient (in m−2 per coccosphere) resulting from an integral over the size distribution of the cores:

acc*(λ)=0P(μ,σ,a0)Qa(k,rc)πrc2da0    (63)

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:

acc(λ)=Ncc[FscFnc+(1Fnc)]acc*(λ)    (64)

where Fnc is the proportion of naked cores and Fsc is the fraction of naked cores that remain in the surface after their liths have been shed (Fsc = number of surviving naked cores /Ncc). 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):

bbpic(λ)=Ncc[(1Fnc)σbbs(λ)+FncσbbNl(λ)]    (65)

Even though cell densities in bloom conditions for E. huxleyi are somewhat arbitrary, a minimal value of 109 cells m−3 has been suggested (Tyrrell and Merico, 2004), with cell concentration exceeding 1011 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 Ncc = 109 and 2 × 1010 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 Qbbl is only a function of λw/ao, as shown in Figure 4. For a coccolith with semi-major axis length ao, 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 Qbbl = ωbsm(1 − RFrg). As the wavelength gets shorter, parts of the surface will start to appear as rough and Qbbl will increase rapidly to reach a maximum value of Qbbl = ωbrgFrg + ωbsm(1 − Frg) 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).


Figure 4. Backscattering efficiency factor Qbbl vs. the ratio of wavelength in air to the length of the semi-major axis of a coccolith ao for (A) various values of the aspect ratio r0 = ao/bo of the elliptical coccolith, and (B) various values of the gap shape factor rsh = Δsgminsgmax.

The backscattering efficiency for an elliptical coccolith with various aspect ratios r0 = ao/bo 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 Frg of rough area to total area because of the increasingly elliptical shape. The effect of the gap shape factor rsh (defined in Equation 9) is shown in Figure 4B for an elliptical axis ratio r0 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 (rsh ≈ 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 (rsh ≈ 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.

Spectral Backscattering Efficiency for Populations of Elliptical Coccoliths

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.


Figure 5. Average backscattering efficiency for a second order shifted Gamma distribution of liths as a function of mean a0 for a standard deviation of 0.15 (A) and as a function of standard deviation for a mean a0 of 1.4 μm (B).

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 σbbNlbbs for a coccolithophore with a single layer of liths and core radius rc 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.


Figure 6. Ratio of free coccolith to coated coccolithophore backscattering cross-section σbbNlbbs, 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.

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 ao 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 bbp (λ) can be obtained by multiplication of σbbtot (λ) with the initial number concentration of coated cores, Ncc. 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-to-cell 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., 1991, 1993; Gordon 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 (Balch et al., 1993), and (iii) the preferential sinking of coccospheres out of the surface layer due to their larger size, but see Monteiro et al. (2016).


Figure 7. Changes in the total spectral backscattering cross-section σbbtot magnitude (A) and spectral shape (B) as a coccolithophore with a core radius rc of 3.7 μm, a lith semi-major axis a0 of 1.7 μm, and 4 layers of liths sheds its liths.

The spectral shape of bbp(λ) is commonly described by the relationship (Gordon and Morel, 1983):

bbp(λ)=bbp(λ0)(λ/λ0)ξ    (66)

where λo is a reference wavelength and ξ represents the dimensionless spectral slope parameter of bbp(λ). 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 with multispectral bbp(λ) measurements in E. huxelyi cultures who reported a steepening of the spectral slope of bbp(λ) 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 above-water remote sensing reflectance Rrs(λ) (in sr−1) during a lith shedding event of E. huxleyi blooms with initial coccolithophore concentrations ranging between Ncc = 109 and 2 × 1010 cells m−3 is depicted in Figure 8. We assumed a coccolith size distribution with a mean a0 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 Rrs(λ) at the initial stage with all liths attached (Fnc = 0) and Figures 8C,D shows the final stage with all liths shed (Fnc = 1) and no living cores remaining (Fsc = 0). Before the onset of lith shedding, the absorption features of the cores clearly influence Rrs(λ) by reducing reflectance in blue wavelengths (~ 0.4 − 0.47 μm), resulting in peak Rrs(λ) at green wavelengths (~ 0.50 − 0.55 μm). When the liths are freed and no cores are left to absorb light, Rrs(λ) increases by a factor of up to four and the peak in Rrs(λ) shifts toward blue-turquoise wavelengths (~ 0.40 − 0.49 μm), where Rrs(λ) shows little spectral variation. This gives a bright milky turquoise coloring of the ocean, typically associated with E. huxleyi blooms. Our modeled Rrs(λ) spectra are also quantitatively consistent with multispectral ocean color satellite observations of Rrs(λ) in coccolithophore blooms which have been shown to peak at 0.49 μm (Moore et al., 2012).


Figure 8. Variability in the spectral remote sensing reflectance Rrs(λ) (in sr−1) during coccolith shedding by E. huxleyi in blooms with various initial coccolithophore concentrations and for a coccolith size distribution with a mean a0 of 1.7 μm and standard deviation of 0.15 μm and coccospheres with four layers of liths. (A) Initial stage with all liths attached (Fnc = 0), (C) final stage with all liths shed (Fnc = 1) and no living cores remaining (Fsc = 0). Spectra in (B,D) were normalized at a wavelength of 0.55 μm. Note that the lith concentration in panels (C,D) is simply 36 times the cell density.

Next we investigate the influence of changes in the fraction of naked cores (Fnc) or surviving cores (Fsc) on Rrs(λ) during different stages of a bloom with initial E. huxleyi coccolithophore concentration of 109 cells m−3. Figures 9A,B shows how Rrs(λ) changes when the coated cores shed their liths and they die after shedding, i.e., Fsc = 0. Even though there are significant changes in the magnitude of Rrs(λ), 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 Rrs(λ) due to lith shedding will, however, be more pronounced for larger cell concentrations. Figures 9C,D illustrates variability in Rrs(λ) due to changes in the fraction of cores that remain in the surface after lith shedding (Fsc from 0 to 1 when Fnc = 1). Such changes result in much stronger spectral variability in Rrs(λ) by virtue of their absorption features (Figures 9C,D).


Figure 9. Spectral above-water remote sensing reflectance Rrs(λ) in absolute values of sr−1 (A,C) and normalized at 0.55 μm (B,D) for various fractions of naked cores (Fnc = 0 to 1, top row), and various surviving core fractions (Fsc = 0 to 1, bottom row). We used a size distribution of free liths of mean a0 1.7 μm and standard deviation of 0.15 μm, an initial coccolithophore cell concentration of 109 m−3 and four layers of liths.

Lastly, Figure 10 illustrates how changes in the mean size of free liths impact the magnitude and spectral shape of Rrs(λ). The magnitude of Rrs(λ) increases with increasing mean size mostly due to increases in the geometric cross section. The spectral shape of Rrs(λ) 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.


Figure 10. (A) Absolute and (B) normalized spectral remote sensing reflectance Rrs(λ) (in sr−1) for populations of free coccoliths of varying mean size, a0. The standard deviation of the size distribution was 0.15 μm.


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, bbp/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 bbp(λ). The spectral shape of bbp(λ) 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:

λonset=nw4Δsgmax=nw42πa0Pn(arao,brao)N(1+rpg),                                    Qbbl(λonset)=ωbsm    (67)

where N is the number of gaps, rpg is the pillar to gap ratio defined in Equation (9), nw is the refractive index of seawater, and Pn 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, rpg, ar/ao, and br/bo, the lith size a0 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:

Qbbl(λ0)=ωbrgFrg+ωbsm(1Frg)    (68)

where Frg 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 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 bbp(λ) are and Frg. Our results suggest that the mean lith size of the population may be obtained from the spectral shape of bbp(λ) 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 R2 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 bbp(λ) would thus also be morphotype-specific. 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 bbp(λ) change surprisingly little. This is due to the quasi-equality of the backscattering cross-sections of the same amount of liths when free or arranged in a coccosphere. The subtle changes in bbp(λ) 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 (Fnc) appears to be feasible if accurate measurements of the spectral slope of bbp (λ) 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 Rrs(λ) 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 bbp(λ) when coccolithophores shed their liths and naked cores die, changes in spectral Rrs(λ) 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 Rrs(λ) by virtue of their absorption. However, similar effects on Rrs (λ) 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 Rrs(λ) 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 Rrs (λ) 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, Rrs(λ) would peak at violet wavelengths and the magnitude of Rrs(λ) 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 × 1010 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 109 m−3 E. huxleyi cores.


Figure 11. Spectral remote sensing reflectance Rrs(λ) (in sr−1) from the free liths shed by Ncc = 109 E. huxleyi cells m−3 in waters with aph(λ) = 0 and varying CDOM absorption, when (A) all cores have died, Fsc = 0 and when (B) all cores survived, Fsc = 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.

In coccolithophore blooms characterized by high concentrations of PIC (>0.003 mol C m−3) the ocean color remote sensing algorithm of Gordon et al. (2001) is used to derive PIC. Because bands at shorter wavelengths often saturate in these bloom conditions, the algorithm uses Rrs(λ) in three bands in the red and near-infrared (0.670, 0.765, and 0.865 μm) to estimate bb from Rrs at 0.670 μm assuming Rrs = 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 Rrs at 0.765 and 0.865 μm much lower than Rrs at 0.670 μm, thereby providing support for the atmospheric correction approach of Gordon et al. (2001). We suggest that our model may be used to improve on the atmospheric correction approach of Gordon et al. (2001).

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.

Roriπab[2πab+2π(a+b)τl4]=21+(1+r)(τla)    (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 Qbb 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 Qbb 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.

Author Contributions

GN and GF conceptualized the study; GF developed the model. Both authors discussed results and wrote the manuscript.


GN is supported by a European Union's Horizon 2020 Marie Sklodowska-Curie grant (no. 749949).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


We would like to thank Dariusz Stramski and Aimee Neeley for providing the spectral absorption coefficient for an Emiliania huxleyi culture.

Supplementary Material

The Supplementary Material for this article can be found online at:


Aas, E. (1996). Refractive index of phytoplankton derived from its metabolite composition. J. Plankton Res. 18, 2223–2249. doi: 10.1093/plankt/18.12.2223

CrossRef Full Text | Google Scholar

Albert, A., and Mobley, C. (2003). An analytical model for subsurface irradiance and remote sensing reflectance in deep and shallow case-2 waters. Opt. Exp. 11:2873. doi: 10.1364/OE.11.002873

PubMed Abstract | CrossRef Full Text | Google Scholar

Armstrong, R. A., Lee, C., Hedges, J. I., Honjo, S., and Wakeham, S. G. (2001). A new, mechanistic model for organic carbon fluxes in the ocean based on the quantitative association of POC with ballast minerals. Deep Sea Res. Part II Top. Stud. Oceanogr. 49, 219–236. doi: 10.1016/S0967-0645(01)00101-1

CrossRef Full Text | Google Scholar

Babin, M., Stramski, D., Ferrari, G. M., Claustre, H., Bricaud, A., Obolensky, G., et al. (2003). Variations in the light absorption coefficients of phytoplankton, nonalgal particles, and dissolved organic matter in coastal waters around Europe. J. Geophys. Res. 108:3211. doi: 10.1029/2001JC000882

CrossRef Full Text | Google Scholar

Balch, W. M., Drapeau, D. T., Bowler, B. C., Lyczskowski, E., Booth, E. S., and Alley, D. (2011). The contribution of coccolithophores to the optical and inorganic carbon budgets during the southern ocean gas exchange experiment: new evidence in support of the “great calcite belt” hypothesis. J. Geophys. Res. 116:C00F06. doi: 10.1029/2011JC006941

CrossRef Full Text | Google Scholar

Balch, W. M., Gordon, H. R., Bowler, B. C., Drapeau, D. T., and Booth, E. S. (2005). Calcium carbonate measurements in the surface global ocean based on moderate-resolution imaging spectroradiometer data. J. Geophys. Res. 110:C07001. doi: 10.1029/2004JC002560

CrossRef Full Text | Google Scholar

Balch, W. M., Holligan, P. M., Ackleson, S. G., and Voss, K. J. (1991). Biological and optical properties of mesoscale coccolithophore blooms in the Gulf of Maine. Limnol. Oceanogr. 36, 629–643. doi: 10.4319/lo.1991.36.4.0629

CrossRef Full Text | Google Scholar

Balch, W. M., Kilpatrick, K. A., Holligan, P., Harbour, D., and Fernandez, E. (1996). The 1991 coccolithophore bloom in the central North Atlantic. 2. Relating optics to coccolith concentration. Limnol. Oceanogr. 41, 1684–1696. doi: 10.4319/lo.1996.41.8.1684

CrossRef Full Text | Google Scholar

Balch, W. M., Kilpatrick, K., Holligan, P. M., and Cucci, T. (1993). Coccolith production and detachment by Emiliania huxleyi (Prymnesiophyceae). J. Phycol. 29, 566–575. doi: 10.1111/j.0022-3646.1993.00566.x

CrossRef Full Text | Google Scholar

Broecker, W., and Clark, E. (2009). Ratio of coccolith CaCO3 to foraminifera CaCO3 in late Holocene deep sea sediments. Paleoceanography 24:PA3205. doi: 10.1029/2009PA001731

CrossRef Full Text | Google Scholar

Brown, C. W., and Yoder, J. A. (1994). Coccolithophorid blooms in the global ocean. J. Geophys. Res. 99:7467. doi: 10.1029/93JC02156

CrossRef Full Text | Google Scholar

Cauchy, A.-L. (1850). Mémoire sur la Rectification des Courbes et la Quadrature des Surfaces Courbes. Mémoires l'académie des Sci. tome XXII, 167–177. Available online at: (Accessed March 26, 2018).

Cook, S. S., Whittock, L., Wright, S. W., and Hallegraeff, G. M. (2011). Photosynthetic pigment and genetic differences between two Southern Ocean morphotypes of Emiliania huxleyi (haptophyta). J. Phycol. 47, 615–626. doi: 10.1111/j.1529-8817.2011.00992.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Fournier, G., and Neukermans, G. (2017). An analytical model for light backscattering by coccoliths and coccospheres of Emiliania huxleyi. Opt. Express 25:14996. doi: 10.1364/OE.25.014996

PubMed Abstract | CrossRef Full Text | Google Scholar

Gordon, H. R. (2006). Backscattering of light from disklike particles: is fine-scale structure or gross morphology more important? Appl. Opt. 45:7166. doi: 10.1364/AO.45.007166

PubMed Abstract | CrossRef Full Text | Google Scholar

Gordon, H. R., Boynton, G. C., Balch, W. M., Groom, S. B., Harbour, D. S., and Smyth, T. J. (2001). Retrieval of coccolithophore calcite concentration from SeaWiFS Imagery. Geophys. Res. Lett. 28, 1587–1590. doi: 10.1029/2000GL012025

CrossRef Full Text | Google Scholar

Gordon, H. R., and Du, T. (2001). Light scattering by nonspherical particles: application to coccoliths detached from Emiliania huxleyi. Limnol. Oceanogr. 46, 1438–1454. doi: 10.4319/lo.2001.46.6.1438

CrossRef Full Text | Google Scholar

Gordon, H. R., and Morel, A. Y. (1983). Remote Assessment of Ocean Color for Interpretation of Satellite Visible Imagery: A Review. Washington, DC: American Geophysical Union.

Google Scholar

Gordon, H. R., Smyth, T. J., Balch, W. M., Boynton, G. C., and Tarran, G. A. (2009). Light scattering by coccoliths detached from Emiliania huxleyi. Appl. Opt. 48:6059. doi: 10.1364/AO.48.006059

PubMed Abstract | CrossRef Full Text | Google Scholar

Groom, S. B., and Holligan, P. M. (1987). Remote sensing of coccolithophore blooms. Adv. Sp. Res 7, 73–78. doi: 10.1016/0273-1177(87)90166-9

CrossRef Full Text | Google Scholar

Hagino, K., Bendif, E. M., Young, J. R., Kogame, K., Probert, I., Takano, Y., et al. (2011). New evidence for morphological and genetic variation in the cosmopolitan coccolithophore Emiliania huxleyi (prymnesiophyceae) from the COX1b-ATP4 genes. J. Phycol. 47, 1164–1176. doi: 10.1111/j.1529-8817.2011.01053.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hedley, J. (2012). Modelling the optical properties of suspended particulate matter of coral reef environments using the finite difference time domain (FDTD) method. Geo-Mar. Lett. 32, 173–182. doi: 10.1007/s00367-011-0265-8

CrossRef Full Text | Google Scholar

Hoffmann, R., Kirchlechner, C., Langer, G., Wochnik, A. S., Griesshaber, E., Schmahl, W. W., et al. (2015). Insight into Emiliania huxleyi coccospheres by focused ion beam sectioning. Biogeosciences 12, 825–834. doi: 10.5194/bg-12-825-2015

CrossRef Full Text | Google Scholar

Holligan, P. M., Fernández, E., Aiken, J., Balch, W. M., Boyd, P., Burkill, P. H., et al. (1993). A biogeochemical study of the coccolithophore, Emiliania huxleyi, in the North Atlantic. Glob. Biogeochem. Cycles 7, 879–900. doi: 10.1029/93GB01731

CrossRef Full Text | Google Scholar

Hopkins, J., Henson, S. A., Painter, S. C., Tyrrell, T., and Poulton, A. J. (2015). Phenological characteristics of global coccolithophore blooms. Glob. Biogeochem. Cycles 29, 239–253. doi: 10.1002/2014GB004919

CrossRef Full Text | Google Scholar

Iglesias-Rodríguez, M. D., Brown, C. W., Doney, S. C., Kleypas, J., Kolber, D., Kolber, Z., et al. (2002). Representing key phytoplankton functional groups in ocean carbon cycle models: Coccolithophorids. Glob. Biogeochem. Cycles 16, 47-1–47-20. doi: 10.1029/2001GB001454

CrossRef Full Text | Google Scholar

Matrai, P. A., and Keller, M. D. (1993). Dimethylsulfide in a large-scale coccolithophore bloom in the Gulf of Maine. Cont. Shelf Res. 13, 831–843. doi: 10.1016/0278-4343(93)90012-M

CrossRef Full Text | Google Scholar

Monteiro, F. M., Bach, L. T., Brownlee, C., Bown, P., Rickaby, R. E., Poulton, A. J., et al. (2016). Why marine phytoplankton calcify. Sci. Adv. 2:e1501822. doi: 10.1126/sciadv.1501822

PubMed Abstract | CrossRef Full Text | Google Scholar

Moore, T. S., Dowell, M. D., and Franz, B. A. (2012). Detection of coccolithophore blooms in ocean color satellite imagery: a generalized approach for use with multiple sensors. Remote Sens. Environ. 117, 249–263. doi: 10.1016/j.rse.2011.10.001

CrossRef Full Text | Google Scholar

Morel, A., and Bricaud, A. (1981). Theoretical results concerning light absorption in a discrete medium, and application to specific absorption of phytoplankton. Deep Sea Res. Part A. Oceanogr. Res. Pap. 28, 1375–1393. doi: 10.1016/0198-0149(81)90039-X

CrossRef Full Text | Google Scholar

Morel, A., and Bricaud, A. (1986). Inherent optical properties of algal cells, including picoplankton. Theoretical and experimental results. Can. Bull. Fish. Aquat. Sci. 214, 521–559.

Google Scholar

Nayak, A. R., McFarland, M. N., Sullivan, J. M., and Twardowski, M. S. (2018). Evidence for ubiquitous preferential particle orientation in representative oceanic shear flows. Limnol. Oceanogr. 63, 122–143. doi: 10.1002/lno.10618

PubMed Abstract | CrossRef Full Text | Google Scholar

Neeley, A. R., Freeman, S. A., and Harris, L. A. (2015). Multi-method approach to quantify uncertainties in the measurements of light absorption by particles. Opt. Express 23:31043. doi: 10.1364/OE.23.031043

PubMed Abstract | CrossRef Full Text | Google Scholar

Neukermans, G., and Fournier, G. (2018). “Rrs_bbp_Ehux: analytical model for remote sensing reflectance and backscattering coefficient from Emiliania huxleyi blooms,” GitHub. doi: 10.5281/zenodo.1222124

CrossRef Full Text

Neukermans, G., Oziel, L., and Babin, M. (2018). Increased intrusion of warming Atlantic water leads to rapid expansion of temperate phytoplankton in the Arctic. Glob. Chang. Biol. doi: 10.1111/gcb.14075. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

Paasche, E. (2001). A review of the coccolithophorid Emiliania huxleyi (Prymnesiophyceae), with particular reference to growth, coccolith formation, and calcification-photosynthesis interactions. Phycologia 40, 503–529. doi: 10.2216/i0031-8884-40-6-503.1

CrossRef Full Text | Google Scholar

Riebesell, U., Bach, L. T., Bellerby, R. G. J., Monsalve, J. R. B., Boxhammer, T., Czerny, J., et al. (2016). Competitive fitness of a predominant pelagic calcifier impaired by ocean acidification. Nat. Geosci. 10, 19–23. doi: 10.1038/ngeo2854

CrossRef Full Text | Google Scholar

Stramski, D., Bricaud, A., and Morel, A. (2001). Modeling the inherent optical properties of the ocean based on the detailed composition of the planktonic community. Appl. Opt. 40:2929. doi: 10.1364/AO.40.002929

PubMed Abstract | CrossRef Full Text | Google Scholar

Tyrrell, T., and Merico, A. (2004). “Emiliania huxleyi: bloom observations and the conditions that induce them,” in Coccolithophores, eds H. R. Thierstein and J. R. Young (Berlin; Heidelberg: Springer), 75–97.

Google Scholar

van de Hulst, H. C. (1957). Light Scattering by Small Particles, 1st edn. New York, NY: John Wiley and Sons.

Google Scholar

Voss, K. J., Balch, W. M., and Kilpatrick, K. A. (1998). Scattering and attenuation properties of Emiliania huxleyi cells and their detached coccoliths. Limnol. Oceanogr. 43, 870–876. doi: 10.4319/lo.1998.43.5.0870

CrossRef Full Text | Google Scholar

Winter, A., Henderiks, J., Beaufort, L., Rickaby, R. E. M., and Brown, C. W. (2014). Poleward expansion of the coccolithophore Emiliania huxleyi. J. Plankton Res. 36, 316–325. doi: 10.1093/plankt/fbt110

CrossRef Full Text | Google Scholar

Young, J., Geisen, M., Cros, L., Kleijne, A., Sprengel, C., Probert, I., et al. (2003). A guide to extant coccolithophore taxonomy. J. Nannoplankt. Res. Spec. Issue. Available at: [Accessed February 5, 2018].

Google Scholar

Young, J. R., Poulton, A. J., and Tyrrell, T. (2014). Morphology of Emiliania huxleyi coccoliths on the northwestern European shelf – is there an influence of carbonate chemistry? Biogeosciences 11, 4771–4782. doi: 10.5194/bg-11-4771-2014

CrossRef Full Text | Google Scholar

Zhai, P.-W., Hu, Y., Trepte, C. R., Winker, D. M., Josset, D. B., Lucker, P. L., et al. (2013). Inherent optical properties of the coccolithophore: Emiliania huxleyi. Opt. Express 21:17625. doi: 10.1364/OE.21.017625

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: coccolithophores, optical model, coccolith morphology, backscattering coefficient, backscattering efficiency, size distribution, ocean color remote sensing, hyperspectral

Citation: Neukermans G and Fournier G (2018) Optical Modeling of Spectral Backscattering and Remote Sensing Reflectance From Emiliania huxleyi Blooms. Front. Mar. Sci. 5:146. doi: 10.3389/fmars.2018.00146

Received: 10 February 2018; Accepted: 10 April 2018;
Published: 07 May 2018.

Edited by:

Astrid Bracher, Alfred Wegener Institut Helmholtz Zentrum für Polar und Meeresforschung, Germany

Reviewed by:

Guangming Zheng, National Oceanic and Atmospheric Administration (NOAA), United States
Michael Twardowski, Florida Atlantic University, United States
Martin Hieronymi, Helmholtz-Zentrum Geesthacht Zentrum für Material- und Küstenforschung, Germany

Copyright © 2018 Neukermans and Fournier. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Griet Neukermans,