Optical backscattering and linear polarization properties of the colony forming cyanobacterium Microcystis.

Light scattering characteristics of the cyanobacterium Microcystis are investigated with numerical models for sphere aggregates. During summer bloom seasons, Microcystis is prevalent in many inland waters across the globe. Monitoring concentrations with remote sensing techniques requires knowledge of the inherent optical properties (IOPs), especially the backscattering properties of Microcystis cells and colonies in natural settings. In situ measurements in waters dominated by Microcystis blooms have previously detected extremely high backscattering ratios, i.e., bb/b>0.043 at 443 nm [1], the highest to our knowledge in the natural environment. These highbb/bvalues could hold promise as a diagnostic tool in identifying and monitoring Microcystis using optical approaches. However, it has been unclear how this type of optically 'soft' organic particle can generate such highbb/bvalues. In this study, the Multiple Sphere T-matrix (MSTM) model is used to calculate the IOPs of model colonies, including bb/b. Colony sizes in the model ranged from several cells to several hundred and both colony packing density and cell gas vacuole content were varied. Results are compared with model results for equivalent-volume spheres (EVS) and direct in situ measurements. Colony formation was required in the modeling to reproduce the high bb/bconsistent with in situ measurements. The combination of moderate to very dense colony (packing density >30%) and high gas vacuole content in individual cells (volume percentage >20%) was the most favorable condition leading to rapid increases in bb/bwith increasing number of cells Ncell of the colony. Significant linear correlations were observed betweenbb/b and Ncell1/3 for these colonies, wherebb/b increased beyond 0.04 once cell number reached about 1000 cells in the case with the most densely packed cells and highest gas vacuole content. Within commonly observed colony sizes (Ncell <106), colonies with high gas vacuole content exhibited bb/bvalues up to 0.055, consistent with direct measurements from Lake Erie. Polarized scattering was also of interest as a diagnostic tool, particularly with future Earth-orbiting polarimeters being deployed for the NASA Plankton, Aerosols, Cloud, ocean Ecosystem (PACE) mission. The Degree of Linear Polarization (DoLP), expressed by the ratio of two Mueller matrix elements-P12/P11, decreased with increasing colony cell number for Microcystis. Another ratio of two Mueller matrix elementsP22/P11, an index for nonsphericity, also decreased with increasing colony size. In addition to higher relative backscattering, greater colony packing density and larger gas vacuole sizes both led to lower DoLP peak magnitude and lowerP22/P11. An optical opposition feature due to constructive phase interference that was observed previously for cosmic dusts is also present for these modeled colonies, manifested by a narrow intensity peak and negative polarization dip near exact backscattering direction, gradually forming as colony size increases.


Introduction
Toxic blooms of the freshwater cyanobacterium Microcystis are found across the globe and are increasing in severity and duration in many inland systems [1][2][3][4][5]. Their ecological success is partly attributable to a tolerance of high light and their ability to move vertically using buoyancy control from internal gas vacuoles and overall cell mass. Microcystis also is a colonial organism, and colony formations can amplify the ascending or descending rates. Even though individual cells are on the order of a few microns in diameter, colonies can reach very large formations on the scale of millimeters [6,7]. Among natural water constituents, phytoplankton have special importance in regulating primary productivity, serving as a reservoir of biomass and carbon, and, in some cases, causing harmful algal blooms (HABs) that may result in hypoxia and/or toxin production. Phytoplankton pigmentation has enabled detection through absorption and fluorescence measurements [8], and detection of chlorophyll from satellites has been possible since the launch of Coastal Zone Color Scanner (CZCS) in 1978. In situ and remote characterization of phytoplankton using light scattering has received far less attention, albeit with important exceptions (e.g., [9][10][11][12]). One of the central challenges in making progress with scattering as a diagnostic measurement for phytoplankton is the lack of satisfactory models, as even the simplest phytoplankton cells are far too complex to precisely simulate with any existing model. However, an interest here is the degree to which relatively simple models that incorporate first order drivers such as bulk refractive index, approximate cell shape, a unique shell layer, presence or absence of gas vacuoles, and colony forming tendencies, may be sufficient to identify specific phytoplankton types within complex particle fields. The use of remote sensing for detecting cyanobacteria in freshwater systems has been increasing in recent years [13]. There are optical features in the remote sensing reflectance spectrum that indicate the presence of cyanobacteria [14][15][16][17][18]. Due to the presence of gas vacuoles, backscatter properties are enhanced in Microcystis (and other vacuolate genera) [11,19]. Modeling studies of individual cells have previously found Microcystis has a backscatter ratio up to 1%, higher than many other species that lack a hard shell [4]. Field observations of Microcystis blooms however measured much higher backscatter ratios [1]. The mechanism for higher observed values is not fully understood and currently unknown from experimental models. It is important to determine and understand the sources of backscatter enhancement for interpretation of optical and remote sensing data.
The backscattering ratio b b /b is an inherent optical property (IOP) of the bulk particle population of a water mass that has been closely linked to particle composition [20][21][22]. For remote sensing applications, the ocean reflectance is given often related to the ratio of backscattering to absorption (e.g., [23]) and recent work suggests b b /b may have a non-negligible influence on reflectance [24]. The Degree of Linear Polarization (DoLP) measures the ability of the particles to change the linear polarization state of incident light. The DoLP angular signal may contain rich information about the size, shape, and composition or internal structure of particles [25]. Measurements or modeling of DoLP for oceanic particles are however scarce. In optics models dating back to the 1970's, oceanic particles have typically been approximated as homogeneous spheres (Lorenz-Mie) or coated spheres (e.g., [26]). With the development of new light scattering algorithms and parallel processing computer hardware, modeling moderate to large non-spherical particles becomes possible. Scattering models for coated spheres [4,[27][28][29][30], spheroids [28,[31][32][33][34], hexahedron [30,35] and specific non-spherical shapes [36][37][38][39][40][41] have now been used to simulate IOPs of marine particles. Non-spherical particles produce different backscattering and linear polarization angular signals than comparable (e.g. equivalent-volume) homogeneous spheres [31,42], although simulations of the polarized scattering characteristics of phytoplankton are scarce. Also, only non-spherical particle shapes can produce the angular curve shape in the Mueller matrix element P 22 /P 11 observed by the Voss and Fry measurement on ocean waters [43]. There is now a need for algorithms to interpret polarized, water-leaving radiances in terms of particle field characteristics for the upcoming NASA PACE mission, which will have two polarimeters [44]. A better understanding of such light scattering properties for Microcystis may enable the use of simple optical property measurements, both from in-water sensors and derived from remote satellite imagers, for environmental monitoring. Pursuing more accurate interpretations of direct measurements of IOPs [1,[45][46][47], particle sizes, geometries and types [48,49] through realistic particle IOP modeling is a first step toward this goal.
The focus of this work is modeling the backscattering and linear polarization properties of Microcystis, a colony forming phytoplankton often responsible for intense, harmful blooms during bloom seasons in inland waters. Microcystis produces the hepatotoxin microcystin, a threat to animal and human health that has been responsible for periodic suspension of large municipal drinking supply systems [50]. In situ measurements have previously detected very high particulate backscattering ratios (with extreme values exceeding 7%) in waters dominated by Microcystis aeruginosa [1,[51][52][53]. These are the highest reported values to our knowledge in natural waters and could possibly hold promise as a diagnostic tool for identifying Microcystis using optical approaches. Individual Microcystis cells have gas vacuoles. Modeling work and measurements in cultures have shown these gas vacuoles elevate backscattering [4,54]. However, in axenic laboratory cultures, it is difficult to get Microcystis to form massive colonies that are similar to those in natural waters except for certain strains [55]. Scattering simulations from single vacuolated cells do not reproduce b b /b high enough to explain in situ measurements [1,52]. In previous modeling studies, single coated spheres were used to model Microcystis [4]. These cells typically aggregate in dense colonies where cell numbers in a single colony can exceed ∼10 6 cells in natural settings. We thus hypothesize that the colony structure may significantly impact the backscattering and polarized scattering properties of this phytoplankton.

Methods
Microcystis colonies are typically comprised of very densely packed substructures with at least ∼10 3 cells ( Fig. 1(a)), which are then interconnected into much larger colony superstructures that can be several millimeters in their maximum dimension [7]. Figure 1(b) shows part of a commonly observed massive colony. Conceptually, a single colony can be seen as an imaginary volume filled by a large number of small particles (cells) in quasi-random positions and random orientations. According to the definition by Mishchenko in [56], it is a type-1 Discrete Random Medium (DRM). Previous studies on the scattering properties of densely packed DRM mainly deal with materials such as soot aggregates [57,58], snow [59,60] or cosmic dusts [61][62][63][64]. Studies on the multiple scattering and coherent backscattering effects of visible light by certain types of cosmic dusts are relevant to this study. In later sections, similar multiple scattering features that were previously observed for cosmic dust aggregates can be observed in the Mueller matrix elements of Microcystis colonies. In our modeling setup, individual Microcystis cells are represented as three-layered spheres. Particle layer complex refractive indices m are given relative to a water refractive index of 1.33. Incident wavelength is 660 nm in air, the wavelength of many commercial scattering sensors. For each individual cell (Fig. 2), the outermost layer is the mucus membrane (m=1.05), the middle layer is the cell cytoplasm (m=1.035+i0.004), and the center is occupied by a gas vacuole (m=0.75) (Fig. 2(a), 2(b)). Values of refractive indices used here are generic for phytoplankton and within the range of numbers reported in previous studies [4,19,54]. The thickness of the mucus membrane relative to the radius of the entire cell (including mucus) is fixed at 0.14, determined from microscopic observations. The volume fraction of the gas vacuole relative to cell (excluding mucus)ρ gas is varied from 0% to 50%. A volume fraction of 10% to 30% is a common value for buoyant cells near the water surface; 50% is considered heavily vacuolated [4]. During summer bloom periods, values up to 90% have been measured [65]. Note Microcystis gas vacuoles are actually clusters of individual cylindrically shaped vacuoles with nominal base diameter of ∼75 nm and height ∼450 nm. Here the gas vacuole cluster is modeled as a single sphere, a constraint of the light scattering algorithm applied. Individual cells are then packed together to form a colony with a certain density, where ρ agg defines the fractional volume of the colony aggregation occupied by the cells, i.e., 1-ρ agg represents the interstitial volume fraction.    Figure 2(a) contains colonies with individual cell diameter D cell =2µm, packing densities of ρ agg =30% and 74%, and their corresponding equivalent-volume sphere (EVS). Each EVS has the same volume as its colony counterpart. The interior of EVS has the same 3-layer composition (mucus membrane, cytoplasm and gas vacuole). By comparing EVS and colony's backscattering properties, we can see the impact of the colony structure on backscattering. In turn, we can see how EVS compares to the backscattering signals of colonies. N cell is the number of cells in the colony. For each column in Fig. 2(ab), the three types of particles have the same total cell volume V, thus the same equivalent volume diameter D eqv = (3V/4π) 1/3 . For colony (2), D max is the diameter of its circumscribing sphere; for colony (3), D max is the longest diagonal of its circumscribing box. Figure 2(b) is the same except colony D max =5µm. Cells 2µm in size are included as this is consistent with the smallest sizes reported for the Microcystis spp. [7,55,66]. For cases with ρ agg =74%, cells are packed according to the Face Centered Cubic (FCC) lattice for equal spheres [67]. This lattice has the highest possible packing density, defined with respect to the circumscribing box volume. As mentioned, fractional gas vacuole volume ρ gas was varied from 0% to 50% in the experiments. For computational efficiency, i.e., to obtain results for highest N cell possible, cell layers in the cell were concentric and cell sizes and internal configuration for a given colony were kept constant. In later sections, the effects of varying cell size, cylindrical vacuole shapes and random gas vacuole positions are assessed.
The Multiple Sphere T-matrix (MSTM) method [68] is used to compute colony optical properties. In this approach, any number of spheres can be inserted in the computational space as long as no spheres intersect. Random orientation averaging for the entire colony is computed. The MSTM method includes analytic orientation averaging [25,69] that is more accurate and efficient than discrete orientation averaging performed in similar methods such as the Discrete Dipole Approximation (DDA) [70]. For equivalent-volume sphere (EVS), the standard Lorenz-Mie algorithm (Matscat [71,72]) is used since it is much faster than the MSTM method. The MSTM algorithms have been used in numerous previous applications. In particular, it was heavily relied upon in the effort to bridge the phenomenological radiative transfer theory to the rigorous Maxwell's equations in the study of electromagnetic scattering by discrete random medium [56].
Three assumptions were made for the model colonies: 1) gas vacuole is a sphere in individual cells; 2) gas vacuole is at the center of the cell; 3) the cells in colonies have equal sizes. Assumption 1 was made because the MSTM algorithm only accepts spheres. Assumptions 2 and 3 were made purely due to computational speed limits of the MSTM. To assess the impacts of assumptions 1 and 2, the Finite Difference Time Domain (FDTD) method [37] was used to model a special case where the gas vacuoles were modeled as non-centered clusters of gas cylinders. In the FDTD modeling setup, discrete orientation averaging was performed. For three colony structures, 60 different orientations were simulated and averaged (180 in total). For assumption 3, a special colony with varying cell sizes was computed and compared to a comparable equal-cell colony.
For a particle assemblage with randomly positioned particles, after orientation averaging, the normalized Mueller matrix (P ij , i, j = 1, 4) is given by [62] [73]: where is the Stokes vector. Subscripts i and s represent the incident and scattered light fields.C sca is the scattering cross section of the particle assemblage, θ is the scattering angle, and R is the distance to the observation point in the far field outside the particle volume. Output from the MSTM, Matscat and FDTD algorithms contain Mueller matrix elements in normalized form, i.e., P ij /P 11 , i, j ≠ 1, 1. The Volume Scattering Function (VSF) is C sca P 11 and the backscattering ratio b b /b is given by: For unpolarized incident light, the Degree of Linear Polarization (DoLP) is −P 12 /P 11 . Colony b b /b, P 11 and −P 12 /P 11 were compared for various configurations. Relationships between colony N cell , b b /b and DoLP were also assessed. Modeled b b /b and DoLP were also compared with in situ measurements from Lake Erie. Aforementioned special cases for testing the assumptions are presented in the final sections of the paper. Simulation of each colony with the MSTM algorithm was conducted on the FAU KOKO cluster across 1 or 2 nodes, where each node has 20 Intel Xeon cores and 128GB RAM. Simulation of the largest colony in this study required about 3 weeks.

Backscattering ratios
For a baseline, b b /b as a function of cell diameter for EVS is shown in Fig. 3. Three gas vacuole fractions ρ gas are shown, 50%, 20%, and 0%. Measurements in [1] reported b b /b>0.04 for water dominated by Microcystis colonies, so 0.04 is set as a threshold for very high relative backscattering in model assessments. EVS with diameter D>2µm generally show b b /b on the order of 0.01 or less (Fig. 3). Gas vacuole size shows a significant impact on enhancing b b /b. However, for D<2.6µm, cells with ρ gas =20% sometimes show higher b b /b than cells with ρ gas =50%. The specific enhancement around D = 2µm for the ρ gas =20% curve is local constructive interference for this particular cell configuration. To help explain this, Fig. 4 shows P 11 of a 2µm cell with 3 different gas vacuole sizes computed with both MSTM and Matscat. All comparisons between the two models for individual cells of any configuration and any size agreed within 0.01% RMSE. The case with ρ gas =20% (middle panel) shows clear enhancement in the backscatter regions compared to the other two cases. This explains the enhancement at D = 2µm, ρ gas =20% in Fig. 3(b). On the other hand, as cell diameter increases, the difference between the three curves becomes smaller ( Fig. 3(a)). Eventually, the three curves converge near 200µm. The biggest relative difference in b b /b between the three gas vacuole percentages is within D eqv range 0.2µm to 20µm. This is also the size range where we will conduct our numerical experiments comparing EVS with colonies. Overall, for EVS, b b /b increases with increasing cell diameter, although values never approach the 0.04 threshold.  Microcystis cells and most phytoplankton cells have low refractive index contrast relative to water, i.e., m is close to 1, and are thus sometimes referred to as optically 'soft' or 'thin'. A gas vacuole within the cell, however, provides a strong refractive index contrast relative to the cell content. With steep refractive index contrast going from m=1.035 to m=0.75 from the cytoplasm to the gas vacuole, total internal reflection will occur at this boundary, thus increasing the number of reflections and relative backscattering. Although gas vacuoles reduce the overall refractive index of the cell, they are quite efficient in reflecting back incident light due to the steep gradient in m. Figure 5 displays b b /b of colonies as a function of their number of cells N cell ; corresponding EVS are also shown. All colonies show higher b b /b than their EVS counterparts except D cell =5µm, ρ gas =0% colonies in Fig. 5(b). For such aggregates, following the logic that backscattering should be controlled more by N cell and that total scattering, dominated by forward-scatter, should be controlled more by the particle cross sectional area correlating with N cell 2/3 , it was hypothesized that a good linear correlation might be found between b b /b and N cell /N cell 2/3 = N cell 1/3 . For colonies with gas vacuoles, positive linear correlations between b b /b and N cell 1/3 were indeed observed. Corresponding least-square fits for colonies in Fig. 5 are listed in Table 1 Fig. 5.  Table 1.
The coefficients of determination R 2 of these least-squares fits are generally high. Colonies with gas vacuoles generally show monotonic linear increases in b b /b; even though the least squares fits do not include the b b /b result at N cell =1, the fitted curves tend to trace back to the starting point at N cell =1. This was especially true for the ρ agg =30% colonies, likely because of their spherical colony shapes. The trends for EVS showed different behavior. For both D cell = 2µm and 5µm, colonies with ρ gas =50% have b b /b values that cross the 0.04 line or rapidly approach this threshold when N cell is still lower than 1000. Table 2    Following these trends, all colonies with gas vacuoles surpass b b /b=0.04 by N cell =10 5 . Thus, the combination of colony formation of cells with gas vacuoles provides results consistent with high b b /b values observed previously in Microcystis blooms [1,45,52,53]. According to Eqs. (3)a and 4a, ρ gas =50% colonies with several millions cells will basically block incident light (b b /b approaching 1). In Fig. 1(b), dense aggregates appear dark because they have enough cells in them to block all the incident light. For a given colony size, i.e., fixed N cell , higher gas vacuole content ρ gas and higher packing densities ρ agg result in higher b b /b. Clearly, ρ gas has significant impact on b b /b, as has been shown previously with a coated Lorenz Mie model [4]. Figures 6 to 8 show the Mueller matrix elements P 11 , −P 12 /P 11 , P 22 /P 11 of small, medium, and large colonies with D cell = 2µm, 5µm and different colony shapes (different ρ agg ). For D cell =5µm colonies, the Mueller matrix elements of ρ agg =30% and 74% are almost identical. So only ρ agg =30% results are shown in Fig. 8. The changes in P 11 as colonies grow larger fit the variations of b b /b in Fig. 5. For colonies with gas vacuoles, increases in P 11 backscatter, i.e., in the 90 to 180 degrees range, with increasing colony size are evident. At the same time, the relative increases in backscattering are being compensated by decreases in the angular range of ∼3 to ∼60 deg. The trend is not always consistent, as P 11 backscatter magnitudes show very slight increases with increasing colony size for the D cell = 2µm, ρ gas =20% cases, consistent with their almost constant b b /b in Fig. 5(a). For the D cell = 2µm, ρ gas =0% cases (Fig. 6 right panels), P 11 backscatter decreases, consistent with decreasing b b /b with increasing colony size in Fig. 5(a). For DoLP (−P 12 /P 11 ) and the index for non-sphericity (P 22 /P 11 ) of colonies with gas vacuoles, magnitudes generally decrease with increasing colony size. Colonies produce similar P 22 /P 11 shape to that measured by Voss and Fry for ocean waters and phytoplankton culture [43] [74]. Comparing Figs. 6∼8 panels horizontally, we see that the gas vacuole size has a Fig. 6. Mueller matrix elements P 11 , −P 12 /P 11 , P 22 /P 11 plotted for D cell =2µm, ρ agg =30% colonies of three different sizes (N cell =8, 88, 512). Three columns correspond to three gas vacuole sizes (ρ gas =50%,20%, 0%). Second row panels focus on backscatter (scattering angle 90 o to 180 o ). significant impact on the Mueller matrix elements. Increase in vacuole ρ gas from 0% to 50% leads to increases in P 11 backscatter. Increasing ρ gas also leads to lower DoLP and P 22 /P 11 . Gas vacuole size also changes the curve shape of these elements, with the most obvious effect on DoLP. For colonies with the two D cell sizes (2 and 5µm), the DoLP shapes corresponding to the three vacuole sizes are significantly different. Moreover, the influence of ρ gas on DoLP is more evident than colony shape. Oscillations in the Mueller matrix elements are common for all colonies, reflecting the simplifications in shape required in the model. The symmetric structure of individual cells, equal cell sizes, and the perfect lattice structure for box-shaped colonies are all factors that tend to intensify these oscillations in the Mueller matrix elements. Afterall, the colony sizes considered here are not large enough (not enough cells in the colonies) to smooth out the oscillations. Certain solar system objects without atmospheres show two opposition optical features within several degrees to the exact backscattering directions [75]: 1) a narrow brightness peak; 2) a negative polarization branch accompanied by a narrow dip. Subsequent modeling studies showed that cosmic aggregates with wavelength-sized dielectric grains can produce such features in their Mueller matrix elements [61][62][63][64]76]. A short explanation is that multiple scattering between equal-sized grains in the dense aggregate contributes to constructive phase interference near exact backscattering direction. The same opposition features were observed for the modeled Microcystis colonies. In Fig. 6 second row of panels, for colonies with gas vacuoles, as colony size increase from 8 cells to 512 cells, the narrow intensity peak and polarization dip gradually form. Also, a narrow peak in the P 22 /P 11 element develops. The opposition features intensify with increasing gas vacuole size. Colonies without gas vacuoles did not produce similar features except for a negative polarization branch in the −P 12 /P 11 element. The opposition features are observed in D cell =5µm colonies (Fig. 8) although they are a bit distorted by the oscillations. This could be a unique feature for Microcystis colonies and perhaps could be used as a possible detection tool.

Comparisons with in situ measurements
Selected model results and in situ scattering measurements are compared in Figs. 9 and 10. Measurement data was collected at station 27 (S27) at 41.717°N, -83.37°E, in Maumee Bay, western Lake Erie, during cyanobacteria blooms in the summer of 2013 (for details see [1]). Time series measurements were made at a depth of ∼ 0.5 m. Microscopic analysis of samples showed Microcystis colonies with individual cell diameters ∼2µm dominated particle fields, but suspended sediment particles were also present.
DoLP modeling results are reproduced from Fig. 6 and 7, with in situ measurements overlaid. Predicted model b bp /b p values (dashed lines) reached the measured value (∼0.055) at S27 within 10 6 cells. For all colony geometries, the fits between model and measured DoLP get better as the gas vacuole content ρ gas increases to 50% and the colony size N cell increases. In measured DoLP of S27, a negative polarization branch near backscattering direction is observed and fits the modeled negative polarization. The current instrument [46] only measures up to 170 degrees so it averages out subtle features within several degrees to exact backscattering. Figure 11 shows a comparison between colonies with equal cell sizes and lognormally distributed cell sizes. Oscillations in the Mueller matrix elements are smoothed out by having varying cell sizes. The DoLP curves are smoothed out but the overall shape is still largely influenced by the gas vacuole size. Moreover, for colonies with gas vacuoles, the local minimum in DoLP curve in the forward direction (i.e., 0 o to 90 o ) are, to a certain degree, similar to results from air bubbles in water. Figure 12 shows a comparison between a colony with cylindrical gas vacuole clusters, which is closer to the actual vacuole shape for Microcystis, and the analog colony with spherical gas vacuoles. Simulation of the colony with cylindrical vacuole clusters was carried out with the Finite-Difference Time Domain (FDTD) approach. There are 13 cells in each colony. Due to computer processing constraints, this was the largest colony that could be modeled with FDTD. Individual cell diameter is D cell = 1.75µm, mucus membrane thickness of the cell is 0.25µm.   Each gas vacuole cylinder is 75 nm in base diameter and 450 nm in height. The cylindrical vacuoles are stacked together in a honeycomb structure, closely following the morphology from microscopic observations. A single spherical gas vacuole with equivalent volume was used for the cells in the MSTM modeled colony, which had a ρ gas =0.11%. For both colonies, gas vacuoles are placed randomly in each cell. Figure 12(b) shows their Mueller matrix elements and b b /b values for the two different modeling approaches. Results for b b /b of the two colonies are close. The major difference is the cylindrical gas vacuoles smooth out the oscillations in the Mueller matrix elements caused by perfectly spherical gas vacuoles. From Figs. 11 and 12, we see the assumptions of concentric spherical cell layers and equal cell sizes in colonies indeed intensify the unnatural oscillations in modeling results. This simulation confirms that a sphere is at least suitable as a first order approximation of gas vacuoles inside a cell and gives credence to the results here, and earlier studies using this approach.

Conclusion
Within the cell size and colony size range considered in this modeling study, we see it is possible to reach the observed b b /b>0.04 for an optically soft organic particle such as Microcystis in colony formations. The combination of moderately dense (ρ agg ≥30%) and high gas vacuole content in individual cells (ρ gas ≥20%) led to rapidly increasing b b /b with N cell of the colony. These colonies generate higher b b /b than their equivalent-volume sphere counterparts. Significant linear correlations between b b /b and N cell 1/3 were found for these model colonies. Within physically feasible, typically observed colony sizes (N cell <10 6 ), model colonies with large gas vacuoles were found to reproduce the extremely high b b /b values (up to 0.08) observed in Lake Erie. At a fixed colony size, increasing gas vacuole size and increasing colony packing density resulted in higher b b /b. As observed previously [4], gas vacuole size has a significant impact on b b /b. The Degree of Linear Polarization decreased with increasing colony size and the nonsphericity indicator P 22 /P 11 also decreases. Increasing colony packing density and gas vacuole size also led to lower DoLP maxima and P 22 /P 11 curves. Gas vacuole size had significant influence on the shape of the Mueller matrix elements. However, this influence is probably due to the artificial, perfectly concentric spherical gas vacuole shapes in each cell. The optical opposition feature that was previously observed for cosmic dusts is present for model colonies. Narrow intensity peak and narrow negative polarization dip near 180 degree backscattering direction gradually manifest as colony size increases. By comparing against measured DoLP, we observed the single spherical gas vacuole in individual cells produced unwanted structure and oscillations. It was further confirmed through a 13-cell test colony comparison using FDTD modeling that the colony with spherical gas vacuole produces much more oscillatory P 11 and DoLP curves than the colony with clusters of cylindrical vacuole clusters. Another special case tested the assumption of equal cell size in a colony. A lognormally distributed colony produced smoother Mueller matrix elements than a colony with equal cell sizes. However, these tests demonstrated the MSTM modeling reproduced the lower frequency, broad shape patterns in the functions along with accurate b b /b. The rapidly increasing b b /b with colony size (or cell number) is a new observation. According to Eqs. (3) and 4, ρ gas =50% colonies with several millions cells will basically block incident light (b b /b approaching 1). For large cyanobacteria colonies, a combination of wave optics and geometric optics is needed. A potential algorithm was reported in [64] and is the subject of future research.