E ﬀ ect of a Coccolithophore Bloom on the Underwater Light Field and the Albedo of the Water Column

: The goal of this work is to study the inﬂuence of coccolithophore blooms on the underwater light ﬁeld and albedo of the water column. A coccolithophore is a single-celled alga with spherical cells surrounded by disk-shaped calcite plates (coccolites), which produce strong light scattering. Because of that, we can observe coccolithophore blooms on satellite ocean color images. We calculated the angular underwater radiance distributions and their integral parameters by the exact numerical method with the input parameters, corresponding to real conditions observed in the Barents Sea and Black Sea. Using the results of the exact calculations, we estimated, for various situations, the accuracy of the approximating formulas applied to the assessment of the water radiance reﬂectance and the di ﬀ use attenuation coe ﬃ cients and we make recommendations for their application. As a ﬁnding of practical importance, we can note the estimate of the accuracy of the widely used Gordon’s formula for the di ﬀ use attenuation coe ﬃ cient; this formula results in large errors under strong coccolithophore blooms. We also mention the interesting and important results concerning the features of the asymptotic regime under such conditions.


Introduction
This work aims to study the influence of coccolithophore blooms on the underwater light field and albedo of the water column. This problem is connected to a planetary-scale phenomenon, because now coccolithophore blooms (CB) are widely distributed in the Atlantic, Pacific, and Indian oceans, and many seas. They cover areas exceeding 100,000 square kilometers and have a significant impact on essential physical and biogeochemical processes [1]. CBs are the most potent producers of CaCO 3 and affect the exchange of carbon dioxide between the ocean and the atmosphere significantly, as well as the greenhouse effect, and, consequently, global climate change. The long-term flux of coccoliths to the ocean floor contributes to the formation of chalk and limestone rocks. The most common coccolithophore species, Emiliania huxleyi, can produce dimethyl sulfide, which is being released into the atmosphere and contributes to the formation of clouds and, thus, changes in planetary albedo [1].
In our work, we have formulated the following goals: (1) Prepare a set of input parameters and software for numerical calculations of the angular distributions and spectral characteristics of underwater radiance and the integral characteristics of the underwater light field depending on the intensity of coccolithophorid blooms and other determining factors.
(2) Perform calculations for the Black Sea and Barents Sea and analyze the results obtained.
(3) Assess the errors of the available approximate formulas for calculating the spectral values of the water radiance and the diffuse attenuation coefficients of the water column by comparing them with the data of numerical calculations, and to present recommendations for choosing the optimal formula depending on the determining factors.
(4) Formulate the main conclusions and to evaluate the scientific and practical significance of the results.
Essential results in most of the directions of our research have been obtained before, but the methods have developed, allowing us to move to our target. Purposeful and meaningful studies of CBs have become regular with the advent of satellite color scanners. Such scanners can detect the coccolithophore blooms by observation of the visible spectrum where coccolithophorids manifest themselves, owing to their strong weakly selective light scattering. In our article, we consider CBs in the Barents Sea and Black Sea. CBs occur regularly in the Barents Sea in July-September and in the Black Sea in May-July. We have been collecting satellite data for these regions since 1998 and conducted regular field research [2][3][4][5][6][7][8].
Satellite observations give us a possibility to evaluate the CB phenological characteristics and to study their spatial and temporal variability. But they do not provide us with full information about the properties of individual cells; We can obtain such data only from direct laboratory determination. We know now that Emiliania huxleyi causes CBs in both the Barents Sea and Black Sea; coccolithophores (also called coccolithophorids) are single-celled algae with spherical cells (15-100 microns in diameter). Strong coccolithophore light scattering is produced by disk-shaped plates (coccoliths) surrounding each cell and consisting of calcium carbonate CaCO 3 . We have not enough information yet about coccolithophore optical characteristics and consider that in Sections 2 and 4.
The satellite ocean color sensors can directly observe only the sub-surface layer in which the seawater optical characteristics depend on the depth. However, some problems need reliable and detailed information about optical conditions in deeper waters. These problems include the condition assessment of the underwater visibility, the exploitation of the systems for underwater observation, their construction and development, and the assessment of light regime on different depths for primary photosynthetic production. The most full and detailed information to satisfy the above needs contains in the underwater light field data. Strictly speaking, the term "light field" is determined as a complete description of the angular radiance distribution at a given point. Along with that, we also consider in our paper the derived quantities, such as the underwater irradiance, radiance reflectance, and diffuse reflectance.
The first measurements of the angular distributions of underwater radiance began before the Second World War [9]. Still, the extensive and complete field studies of the underwater light fields with acceptable accuracy became possible only with the advent of modern electronic photodetectors (see, as an example, the electro-optical camera system [10]). In the late 1960s, V. Timofeeva began to conduct model studies of underwater light fields in a specially equipped experimental pool, using milk media with the addition of dye and sorting the concentrations out to provide the required optical characteristics [11].
In 1936, the Russian physicist A. A. Gershun published the first fundamental work on light field theory [12]. He first introduced a concept of the light vector and derived an equation linking its divergence in a turbid medium with light absorption; now, this equation is widely used and known as a Gershun equation. Gershun also constructed a simple underwater radiance meter, currently named after him, the Gershun tube photometer [13].
Since the 1960s, both theoretical and experimental studies of the underwater light fields have developed rapidly; eleven monographs related to the above problem have been published [9,[12][13][14][15][16][17][18][19][20][21]. Comparing the theoretical and experimental methods, we can note some advantages of the former: an ability to get more detailed and accurate results and to perform the calculations in cases when it is challenging to apply the experimental methods; the former is usually cheaper than the latter.
However, to realize these advantages, a theoretician must first have reliable information about the seawater's optical properties, its vertical distribution, and the observation conditions. Secondly, the calculation method must be validated and appropriate to the problem. The researcher must have a sufficiently powerful computer to perform the required calculations. Currently, in most cases, these requirements are provided; in the case of modeling, the first requirement is set by the researcher himself or the customer, depending on the task.
The calculation of the underwater light field is based on the solution of the radiation transfer equation (RTE). This integrodifferential equation has no analytic solution, but effective numerical methods now exist. A description of strict and various approximate methods for RTE solving is given in [19].
We consider three of the most widely used numerical methods: Monte Carlo, Discrete Ordinates Radiative Transfer (DISORT), and Hydrolight. Monte Carlo (M-C) was first used to calculate the propagation of solar radiation in the atmosphere-ocean system [22]. This method applies to any experimental geometry and the spatial distribution of optical characteristics, but it requires much machine time for calculation.
The improved method of discrete ordinates, DISORT, implements RTE-solving in a plane-parallel environment, such as the atmosphere. The works [23] on a smooth interface and [24] on the rough sea surface generalize DISORT for the ocean-atmosphere system.
Hydrolight is a commercial program for numerical RTE solutions in a plane-parallel medium [25]. It allows for calculating the spectral radiance distribution at a given depth in the water column, both for downward and upward radiation fluxes, taking into account the reflection from the bottom and under and above the sea surface. The book [20] presented details of the method used for RTE. The advantage of the Hydrolight software is the high processing speed, saving the same accuracy, as compared with M-C. Due to the substantial elongation of the seawater phase function in the forward direction, to ensure the necessary skill for backscattering, the calculation requires a considerable number of photons (10 7 or more). A simplified version of the Hydrolight, called EcoLight, calculates only the irradiance and does it 30-100 times faster than Hydrolight [25].
The authors of [26,27] used Hydrolight to estimate errors in the measurements of the spectral absorption and attenuation coefficients by AC-9 instruments under conditions of intense coccolithophore blooms (the value of the scattering coefficient b reached 4 m −1 ). The values of the absorption and scattering coefficients a and b, measured under such conditions, were inconsistent with the measured data on the spectral values of the radiance coefficient R = L u /E d . The calculations by the Hydrolight method allowed the authors to estimate the errors that occurred.
In [28], the authors presented the results of the Monte Carlo modeling, providing precise calculations of the water radiance reflectance just beneath the sea surface and the ocean albedo during coccolithophore blooms. As CB effects, they noted more pronounced stratification, the net cooling of the water column, and decreased total water column productivity.
The structure of this article corresponds to the presented goals. This paper has five sections: 1. Introduction; 2. Materials and methods; 3-4. Results obtained; 5. Discussion and Conclusion. Sections 3 and 4 present two different aspects of our work: 3. Effects of the coccolithophore bloom on the underwater light field (3.1), the water radiance reflectance (3.2), and the diffuse attenuation coefficient (3.3); 4. The influence of the coccolithophore bloom on the accuracy of the approximate formulas. Section 5 includes the formulated conclusions on the work and the significance of the results. The list of references counts 65 references.

Field and Satellite Data for Modeling
In our work, we did not set out to reproduce the results of field measurements through numerical modeling; we used these results to approximate the input modeling parameters to the real conditions observed in the Barents Sea and Black Sea. The Shirshov Institute of Oceanology (SIORAS) conducts regular complex expeditions in these seas; one of the main tasks is to study coccolithophore blooms based on satellite and field measurements. In the Barents Sea, such expeditions have been conducted  The results of the research in the Baltic, Norwegian, and Barents Seas in the 2014-2016 cruises demonstrate the effectiveness of an integrated approach to the use of optical methods combining the accuracy of contact ship measurements and the full coverage of the studied area with satellite images [29].
The 2017 voyage, the longest (51 days), proved to be the most fruitful; the measurements were carried out both at the drift stations and continuously on the ship's course using a flow system [7]. At the drift stations, a submerged transparency meter measured the vertical profiles of seawater temperature and the beam attenuation coefficient; a laser spectrometer in the ship laboratory, measured the fluorescence spectra of colored dissolved organic matter (CDOM) and phytoplankton pigments in the water samples, taken from selected horizons; an Integrated Cavity Absorption Meter (ICAM) with an integrating sphere measured the absorption spectra of seawater, filtrates, and suspended particles [30,31]. A flow-through fluorometer with appropriate sensors, installed in a flow system, provided continuous measurements of the fluorescence of chlorophyll a (Chl a), CDOM, sea surface temperature (SST), and salinity (SSS). From the light measurements during daylight hours, we obtained data on the spectral upwelling radiance just beneath the sea surface L u , spectral surface irradiance E s , the underwater Photosynthetically Available Radiation (PAR) values by a Li-COR device, and spectral values of the surface and underwater irradiance with E-meter and RAMSES devices (see Figure 1).
We took the data obtained at drift station 5580 north of the Kola peninsula (70.14 • N, 35.27 • E), the cruise AMK-68 2017 as the reference values for the Barents Sea input modeling parameters. As seen from the above text, we did not measure all parameters required for the numerical calculations. Still, we measured enough of them to derive the rest by using the appropriate models. We describe the algorithms used for that in Section 2.3. There we present the measured data which allowed the modeling mentioned above.
A spectral absorption meter ICAM provided the spectral absorption coefficients of the seawater (from the samples taken from different depths), a(λ); the filtrate (after filtration through a nuclear filter of 0.4µ), a f (λ), and the particulate matter, a p (λ) (determined as the difference between the spectral values of the absorption coefficients of the seawater and the filtrate); a p (λ) = a(λ)−a f (λ). The difference between the values of the filtrate and pure water absorption coefficients allowed for the estimation of the yellow substance (CDOM) absorption coefficient a g (λ) = a f (λ)−a w (λ).
A submerged transparence meter PUM-A measured the vertical profiles of the beam attenuation coefficient c(530) at 530 nm with error 0.005 m −1 down to 200 m. After installing the plastic cuvette in the measuring channel, PUM-A can measure seawater samples in laboratory mode, in particular, in real-time, if it is involved in the flow system. Having the values of the attenuation and absorption coefficients, we obtained the scattering coefficient at 530 nm: b(530) = c(530)−a(530).
The measured data of the apparent optical characteristics, such as the spectral radiance reflectance ρ(λ), and the diffuse attenuation coefficient K d (λ), allowed us to calculate the values of the inherent optical characteristics.
As seen from Figure 1, we had three devices to measure the downward and upward underwater irradiance that are useful for increasing the completeness of the information obtained and its accuracy. This is very important for the validation of the results of the numerical calculations.
A spectral underwater irradiation meter RAMSES ACC-VIS (TRIOS, Germany) has 190 spectral channels in the spectral range of 320-950 nm (spectral resolution is 3 nm) and an accuracy of ±5%. To conduct measurements in the sounding mode, we used a special frame and the data was transmitted to the ship laboratory by a cable. TRIOS provided calibration of the instrument to the National Institute of Standard and Technology (NIST) standards.
The other irradiance meter, E-meter (BIC 2100), has two measuring modules: the deck and underwater, which simultaneously measure the irradiance in four spectral bands 443, 490, 555, and 625 nm. The specialists from SIORAS modernized the device to perform measurements in an autonomous regime, without a cable. We calibrated the instrument at the SIORAS before each expedition, according to the working standard of the spectral irradiance density.
The PAR meter Li-192SA (Li-COR Company, USA-Canada) measured the surface and underwater quantum irradiance with an accuracy of ±5%. The underwater module included two submerged Li-192 photodiode sensors for measuring the downward and upward irradiance. The Li-190SA deck sensor measured the surface PAR. Our modernization made it possible to work autonomously without cables and with continuous recording of depth.
During our sea expeditions, we collected satellite data when the weather allowed us to get them. We used data from different satellite ocean color sensors, such as SeaWiFS, MODIS, OLCI, VIIRS, and AVHRR. Still, for constructing the long series data, the core information was from SeaWiFS and MODIS-Aqua, and we are now going to focus on the OLCI data [7].

Calculation Method
We calculated the angular and spectral radiance distributions, the downward and upward irradiances, and others for real atmospheric and observation conditions, solving the radiative transfer equations by the matrix operator method developed in [32]. In our work, we use a modification [24] of this method.
We consider the ocean-atmosphere system as a set of plane-parallel homogeneous layers. To solve RTE (without taking into account the effects of polarization), we should specify the absorption coefficient of sea water a(λ) and the volume scattering function (VSF) β(θ,λ), where θ is the scattering angle and λ is the wavelength.
For the atmosphere, we use a two-layer model: the upper layer with Rayleigh optical thickness [33] and Rayleigh phase function, and the lower layer with the model phase scattering function by Gordon and Castaño [34] and the optical thickness τ = τ a (869) 869 λ A , where Angstrom index A and the aerosol optical thickness τ a (869) are taken from the data of satellite color scanners.
The surface is considered as a separate layer. Formulas for a smooth surface are presented below; for a wind-roughened surface, we used the results of [35,36].
The layers are numbered downward from 1 to N; their optical properties are completely determined by optical thickness ∆τ i = c i· ∆z i , single scattering albedo ω 0ι , and the phase function p i (θ). The radiances and irradiances can be calculated at the levels numbered as 0, 1, . . . N, the i-th level being a boundary between the layers i and i+1. Hence, level 0 corresponds to the top of the atmosphere (TOA), and level N is the bottom. The version of the matrix operator method used here is based on recurrence formulas (see below) applied to single-layer reflection and transmission operators. A distinctive feature of the algorithm is the ability to use the results obtained by different numerical methods in one calculation. The result for a multilayer system that includes (or does not include) a surface can be assembled from pre-calculated elements.
For each layer we define two reflectances, r i (Ω, Ω 0 ) and r i (Ω, Ω 0 ), and two transmittances t i (Ω, Ω 0 ) and t i (Ω, Ω 0 ). By definition, r i (Ω, Ω 0 ) is the intensity of upwelling radiation in the direction Ω = (µϕ) if the incident downwelling radiance is L(Ω) = δ(Ω, Ω 0 ) (µ = cos θ; θ and φ are the zenith and azimuthal angles of the incident light, µ is always positive here). The definition of transmittance is obvious; a tilde corresponds to the case of upwelling incident radiation in the water and atmosphere layers r i (Ω, Ω i ) = r i (Ω, Ω i ) and t i (Ω, Ω i ) = t(Ω, Ω i ) This does not apply to the surface; for example, for a smooth surface we have the formulas where r fr (µ,n w ) is the Fresnel reflection coefficient, n w is the water refractive index, and t fr (µ,n w ) = 1− r fr (n w ).
From elementary probabilistic considerations, it follows that the formulas for the upwelling radiances U n and the downwelling radiances D n at level n can be expressed in terms of the transmittance T n and the reflectance R n for truncated systems. In operator form, where T n is the transmittance of the system including layers 1,2, . . . n, R n is backward reflectance for the same system, R n -reflectance for layers n+1, n+2, . . . N, R N is the bidirectional reflectance distribution function of the bottom, and1 is the unity operator δ(µ−µ 0 )δ(ϕ−ϕ 0 ). To calculate the truncated operators, we use the following recurrence formulas: Note that the transmittance operator T n is not involved in calculation of the top of the atmosphere (TOA) radiance L t ≡ R 0 .
Equations (2)-(5) are a complete system of recurrence formulas for calculating the required radiances. The initial condition to Equation (2) is the bottom bi-directional reflectivity R N , the initial conditions of the coupled system of Equations (3) and (4) are T 0 =1 and R 0 = 0. The radiance at the TOA can be calculated in a single pass of N steps, starting from the bottom and ending with the TOA. To find the radiances at other levels, a second pass from the top layer down is required.
To obtain a working algorithm, the direct radiance should be separated from the diffuse one. Let us denote r n = r dir n + r n , t n = t dir n + t n , R n = R dir n + R n , T n = T dir n + T n (6) where the operators denoted by the superscript dir are singular, i.e., proportional to a delta function, and boldface characters denote regular (not containing any delta functions) parts of the operators. Substituting this into Equations (3)-(5), one can obtain a system of Fredholm integral equations of the second kind with a kernel yielding a compact operator. Such equations can be easily solved using a standard numerical technique. Assuming that all of the above operators are axially symmetric, we can expand them into a Fourier series in the parameter ϕ-ϕ 0 , and obtain a decoupled system of recurrence formulas of the same form, with integration not over Ω but over the one-dimensional parameter µ.
To find the solution, it is sufficient to choose a discretization method, i.e., a set of nodes and weights in a quadrature formula. Then the problem reduces to solving a system of linear algebraic equations.
For atmospheric layers, we use the Gauss quadrature; for water layers, a set of nodes and weights proposed in [23]; For the surface, a combination of both.
To calculate the one-layer reflectance and transmittance, two methods were used: the discrete ordinate method [37] or Monte Carlo method. The discrete ordinate method and the program DISORT are described in detail in a series of publications by K. Stamnes et al. [38,39]. We use the part of the code intended for the one-layer system [38], simplified and modified so that it calculates and stores full matrices, rather than radiances for a limited set of user-defined angles. The result of the calculation can be stored on a hard disk and used afterwards as an imported one-layer operator. The Monte Carlo method has been included mainly for checking purpose (see Section 5.1).

Input Parameters for the Modeling
As pointed out in Section 2.2, for the numerical calculation of the underwater radiance, we need a knowledge of the seawater absorption coefficient a(λ) and the seawater volume scattering function (VSF) β(θ,λ) (where θ is the scattering angle and λ is the wavelength). At that, we have to know the above parameters not only for the seawater as a whole, but for each of its optically active components. To bring our model to reality, as close as possible, we consider not only the plated cells and coccoliths, but the non-coccolithophore components, such as particles of another origin, colored organic matter, and of course the water itself. We used the parameters, based on the real optical and biological data, measured in our sea expeditions. Table 1 shows the origin of the user settings. Below we present the results obtained at drift stations in the Black Sea and Barents Sea, a procedure for the calculation of the coccolithophore optical characteristics is described afterwards.  Backscattering coefficient-b b (λ) from ρ(λ); b(λ), β(λ), P(λ)-using model [40][41][42][43] Coccolithophore Parameters plated cell and coccolith concentration N coc and N lith -from direct determination or by setting for model; Optical characteristics-from [4,44].
For the Black Sea, we used the data from two stations, both in the northeastern part of the Black Sea near Gelendzhik,~45 • N, 38 • E, at a depth of 30 m [31]. See Station 1_08.06.2017 as an example of a very intensive coccolithophore bloom with a concentration of coccolithophores at N coc = 8.2 × 10 6 cell/L, and Station 1_08.06.2018 with no CB for assessing the effect of the "non-coccolithophore" particles (see below). At Station 1_08.06.2017, the favorable weather conditions (the solar zenith θ 0 was 25 • , a calm sea) supported both our field measurements and the satellite observations. From MODIS satellite data, the aerosol optical thickness τ a (869) equaled 0.146 and the Angstrom parameter 1.70.
We assumed the total absorption coefficient was a sum of the contributions from pure water and the yellow substance (CDOM), neglecting the phytoplankton absorption, because, according to the ICAM data (see Section 2.1), it made an insignificant contribution. As the absorption coefficient for pure water, we took the data by [45] and the spectral CDOM absorption slope S yel = 0.017 nm −1 for the range 400-500 nm and 0.011 nm −1 for 500-700 nm [43].
As the absorption input parameter, we used the value of the seawater absorption coefficient a(440) = 0.09 m −1 at 440 nm, according to the ICAM data [31]. For studying the absorption effect on the light field and its integral parameters, we performed our modeling not only with the above input parameter but also with much higher inputs, up to 0.53 m −1 .
As the input parameter of the coccolithophore bloom, we used the plated cell concentration N coc , and varied it from 0 to 10 7 cell/L. For the ratio of the coccolith concentration N lith to the plated cell concentration N coc (N lith /N coc ), we assumed it to be equal to 35, according to measurements at the station 1_08.06.2017.
We also took into account the existence of the "non-coccolithophore" particles supplied by river runoff [4]; their optical parameters were calculated by using a model [40][41][42][43]. This model was built on the basis of the directly measured VSF data and validated by the results of the analysis of marine particles by microscopic and chemical studies and the retrieval of the size distribution of marine particles from the VSF. Using the model, we estimated the volume concentrations of small (radiuses less than 1 µ) and large (≥1 µ) "non-coccolithophore" particles from the measured data of 2018 under a no coccolithophore bloom condition [31], and obtained v s = 0.14 cm 3 m −3 for the fine fraction and v l = 0.78 cm 3 m −3 for the coarse fraction. The concentration of "non-coccolithophore" particles was assumed to be independent from the coccolithophore concentration and constant.
For the Barents Sea, we derived the model parameters based on the data of 15.07.2017 at station 5580 of the 68th cruise of R/V "Akademik Mstislav Keldysh". The measured concentration of the coccolithophore plated cells was approximately N coc = 5 × 10 6 cell/L, the detached coccoliths N lith = 1.5 × 10 8 cell/L.
After calculating the b coc , b lith , b bcoc , and b blith , we obtained the b n and b bn values of the non-coccolithophore particles: where the b and b b values are from the measured data (see Table 1). The obtained b n and b bn values allowed us to find the parameters v s and v l of the two-parameter model of the scattering properties [42,46]. Performing the calculation, we got v s = 0.35 cm 3 m −3 and v l = 4.62 cm 3 m −3 .
As to the calculation of the coccolithophore optical characteristics, we can mention several publications which considered the results of theoretical calculations and experimental studies of the optical characteristics of coccoliths and plated cells [47][48][49][50][51][52]; we discuss them in Section 5.2.
To create a set of the input parameters for a calculation of the coccolithophore optical characteristics, we used the data of the work [51] which presents an almost full set of the optical characteristics, except polarization properties, obtained from comprehensive laboratory measurements. These measurements provided spectral data on the beam attenuation cross sections (m 2 particle −l ) for coccoliths, plated cells, and naked cells, and the values of the volume scattering function (VSF) at six wavelengths. They measured VSF at every degree between 10 • and 170 • and, for calculating the b b values, extrapolated the measured data from 170 • to 180 • by assuming VSF(θ) for 170-180 • was constant. They estimated the error in the calculation of b b as ± 1%.
For our modeling, we need the whole phase scattering function (PSF) from 0 • to 180 • and to get it, we analyzed the data of [49] which comprise the Mueller scattering matrix, backscattering probability, and depolarization ratio for Emiliania huxleyi. For their numerical calculations, the authors used a realistic non-spherical model, based on electron micrographs of coccolithophore cells, and the discrete dipole approximation. They presented the phase scattering functions for different cell sizes with the equivalent diameters of the cell from 3.75 to 6.75 µ for the plated cells and from 1.63 to 4.73 µ for coccoliths. We compared their calculated functions and different combinations of them with the VSF from [51], sorting out the appropriate particle concentration. Figure 2 shows the final result as compared with the phase scattering functions for the Black Sea and Barents Sea non-coccolithophore particulates. particles by microscopic and chemical studies and the retrieval of the size distribution of marine particles from the VSF. Using the model, we estimated the volume concentrations of small (radiuses less than 1 μ) and large (≥1 μ) "non-coccolithophore" particles from the measured data of 2018 under a no coccolithophore bloom condition [31], and obtained vs = 0.14 cm 3 m −3 for the fine fraction and vl = 0.78 cm 3 m −3 for the coarse fraction. The concentration of "non-coccolithophore" particles was assumed to be independent from the coccolithophore concentration and constant.
For the Barents Sea, we derived the model parameters based on the data of 15.07.2017 at station 5580 of the 68th cruise of R/V "Akademik Mstislav Keldysh." The measured concentration of the coccolithophore plated cells was approximately Ncoc = 5 × 10 6 cell/L, the detached coccoliths Nlith = 1.5 10 8 cell/L.
After calculating the bcoc, blith, bbcoc, and bblith, we obtained the bn and bbn values of the noncoccolithophore particles: bn = b -(bcoc + blith), bbn = bb -(bbcoc + bblith), where the b and bb values are from the measured data (see Table 1). The obtained bn and bbn values allowed us to find the parameters vs and vl of the two-parameter model of the scattering properties [42,46]. Performing the calculation, we got vs = 0.35 cm 3 m −3 and vl = 4.62 cm 3 m −3 .
As to the calculation of the coccolithophore optical characteristics, we can mention several publications which considered the results of theoretical calculations and experimental studies of the optical characteristics of coccoliths and plated cells [47][48][49][50][51][52]; we discuss them in Section 5.2.
To create a set of the input parameters for a calculation of the coccolithophore optical characteristics, we used the data of the work [51] which presents an almost full set of the optical characteristics, except polarization properties, obtained from comprehensive laboratory measurements. These measurements provided spectral data on the beam attenuation cross sections (m 2 particle −l ) for coccoliths, plated cells, and naked cells, and the values of the volume scattering function (VSF) at six wavelengths. They measured VSF at every degree between 10° and 170° and, for calculating the bb values, extrapolated the measured data from 170° to 180° by assuming VSF(θ) for 170°-180° was constant. They estimated the error in the calculation of bb as ± 1%.
For our modeling, we need the whole phase scattering function (PSF) from 0° to 180° and to get it, we analyzed the data of [49] which comprise the Mueller scattering matrix, backscattering probability, and depolarization ratio for Emiliania huxleyi. For their numerical calculations, the authors used a realistic non-spherical model, based on electron micrographs of coccolithophore cells, and the discrete dipole approximation. They presented the phase scattering functions for different cell sizes with the equivalent diameters of the cell from 3.75 to 6.75 μ for the plated cells and from 1.63 to 4.73 μ for coccoliths. We compared their calculated functions and different combinations of them with the VSF from [51], sorting out the appropriate particle concentration. Figure 2 shows the final result as compared with the phase scattering functions for the Black Sea and Barents Sea noncoccolithophore particulates. As seen from Figure 2A, the phase scattering functions for the plated cells and coccoliths are close to each other, as noted in [51], and they are not as elongated as the functions for non-coccolithophore particulates in Figure 2B. The comparison between the integral parameters for the coccolithophore systems (including both the plated cells and coccoliths) and for the non-coccolithophore particles in the Black Sea confirms the observed difference: the backscattering probability b b = b b /b at 555 nm equals 0.027 for the former and 0.017 for the latter; the values of the average cosine of the single scattering are equal to 0.892 and 0.922, respectively. The most elongated PSF for the in non-coccolithophore particles in the Barents Sea was b b (555) = 0.0091, < cos θ ≥ 0.949. Figure 3 presents the spectral values of the absorption and the particle backscattering coefficients for the Black Sea and Barents Sea used in our numerical simulations.   There is a significant difference between the Black Sea and Barents Sea in their values of the absorption coefficient at a wavelength less than 550 nm-the values at 440 nm equal 0.09 and 018 m −1 , respectively (see Table 2 below). A similar difference is also seen in Figure 3C,D between the values of the backscattering coefficient for N coc = 0: 0.008 and 0.020 m −1 , respectively ( Table 2).
The set of the optical characteristics in Table 2 looks nice for the modeling purposes because it presents a wide variety of the ratio between two the main processes determining the light propagation in the water medium: the scattering and absorption; the ratio b/a at 440 nm changes by a factor of about 70 at 555 nm-more than a factor of 22.  (555), the total scattering and backscattering coefficients at 555 nm, m −1 ; ω 0 (440) and ω 0 (555), the single scattering albedo at 440 and 550 nm (non-dimensional); b b (555), the backscattering probability; <cos>, the average cosine of the single scattering (non-dimensional)

Effect of the Coccolithophore Bloom on the Underwater Light Field
The features of the underwater light fields depend on the seawater optical characteristics and their vertical distribution and the illumination conditions (solar zenith and azimuthal angles, atmosphere condition, cloudiness, sea surface state). We considered the general topics of the problem, including the approaches, in previous sections, now we present the results of our numerical calculations for different conditions and discuss them. Figures 4-8 show, for the Black Sea and the Barents Sea, the angular distributions of the underwater radiance at the plane of the sun, depending on the coccolithophore concentration N coc and the absorption coefficient a(440) for different observation conditions, such as the solar zenith angle θ 0 , the cloudiness, and the wind speed.
As already noted in the introduction, the research of laws and features of underwater light fields has been going on for several decades, including field and laboratory studies (with artificial media, imitating real seawater optical properties) and modeling. The last has several important advantages (see Introduction).
To date, there are general ideas developed about solar radiation propagation in the water column. One of these is the occurrence of the asymptotic radiance distribution in deep waters. Figures 4-8 display this feature, and later we will pay special attention to it. We start with the subsurface horizon, just beneath the sea surface (commonly referred 0 − , as opposed to 0 + , just above the sea surface). The red curves in Figures 4-8 present the distributions at horizon 0 − . We see a classic picture showing the so-called Snellius circle, which corresponds to radiation entering the water column through the sea surface. This radiation propagates inside a cone with a cone angle of approximately 48 • , corresponding to the maximum refraction angle in the absence of surface waves (wind speed is 0). We see the Snellius circle better in the bloom absence (N coc = 0); a strong scattering washes it out even directly below the surface.  sunlight-the maximum of the angular distribution corresponds to the solar zenith angle after refraction when crossing the sea surface. At large depths, the angular distribution gradually becomes symmetrical about the vertical axis and unchangeable with further increasing depth. Such a distribution depends only on the seawater optical characteristics at these depths and not at all on the illumination conditions. This is a so-called asymptotic regime and we discuss it more detail in Section 5.3.  distributions in Figure 4 (left) and Figure 6 (left) between the sun angles of 25° and 60° for the Black Sea shows a significant difference only for subsurface horizons and in the absence of blooming. Figure  6 (on the right) shows that cloud conditions significantly affect both the shape of the angular distribution and the absolute radiance values. The changes in the sea surface state due to an increase in wind speed has a significant effect on the crossing of the sea surface by direct sunlight and near the maximum (Figure 7, θ = 25°).  Table 3, we can note the effect both of the scattering and absorption. The absorption causes a vertical elongation of the angular distributions, increasing it (we characterize this parameter as the ratio of the normalized radiances for 0° and 90°). This effect is more significant for the Black noticeably higher at 410 nm than at 670 nm for Ncoc = 10 7 cell/L. We can explain the observed changes by the behavior of the ratio of the backscattering and absorption coefficients bb/a that is seen in Figure  3.     At small depths, we observe the asymmetry of distributions associated with direct sunlight-the maximum of the angular distribution corresponds to the solar zenith angle after refraction when crossing the sea surface. At large depths, the angular distribution gradually becomes symmetrical about the vertical axis and unchangeable with further increasing depth. Such a distribution depends only on the seawater optical characteristics at these depths and not at all on the illumination conditions. This is a so-called asymptotic regime and we discuss it more detail in Section 5.3.
Comparing Figures 4 and 5, one can see the similarity of the presented results, despite the notable difference in the optical characteristics of the considered seas and the solar zenith angles. Figure 5 shows the analogical angular distributions of the underwater radiance for the Barents Sea. The difference between distributions is more pronounced in the absence of blooming (N coc = 0), it is displayed in the shape of the distributions and the absolute radiance values. The comparison of distributions in Figure 4 (left) and Figure 6 (left) between the sun angles of 25 • and 60 • for the Black Sea shows a significant difference only for subsurface horizons and in the absence of blooming. Figure 6 (on the right) shows that cloud conditions significantly affect both the shape of the angular distribution and the absolute radiance values. The changes in the sea surface state due to an increase in wind speed has a significant effect on the crossing of the sea surface by direct sunlight and near the maximum (Figure 7, θ = 25 • ).
Analyzing the angular radiance distributions and absolute values of the underwater radiance in Figures 4 and 5 and Table 3, we can note the effect both of the scattering and absorption. The absorption causes a vertical elongation of the angular distributions, increasing it (we characterize this parameter as the ratio of the normalized radiances for 0 • and 90 • ). This effect is more significant for the Black Sea, where the initial absorption is less (for the Black Sea, N coc = 5 × 10 6 cell/L, this ratio increased by the factor 1.5, for the Barents Sea, only 1.2). It is the same for the asymptotic diffuse attenuation coefficient k: it increased for a(440) = 0.53 m −1 by a factor of 1.7 for the Black Sea (N coc = 5 × 10 6 cell/L), and only 1.2 for the Barents Sea (see Table 3).  The results of our calculations allow us to get information about the effect of coccolithophore blooms on the satellite ocean color observations. Figures 4-8 comprise the data on the upwelling angular radiance distributions within the angle ranges of −180 • ÷ (−90 • ) and 90 • ÷ 180. The red curves in these figures provide us with the useful information about the impact of different factors on the upwelling radiance: the coccolithophore bloom intensity, the CDOM absorption (Figures 4 and 5), the effects of the solar zenith angle and cloudiness (Figure 6), and the sea surface state determined by wind speed (Figure 7). Figure 8 shows the influence of the coccolithophore bloom intensity and the CDOM absorption on the angular and spectral dependencies of the subsurface radiance reflectance r rs . The observation at nadir corresponds to the viewing angle of ±180 • , and we can note that a CB with N coc = 10 7 cell/L results in increasing r rs by about an order of magnitude at 410 nm, and by about a factor of 20 at 670 nm as compared with no bloom. We can also observe a weak angular dependency of r rs at 410 nm for the case of N coc = 10 7 cell/L. In both cases, N coc = 0 or 10 7 cell/L, the r rs value is much higher at 410 nm than at 670 nm; the reason for that is that there is a significant difference between the values of the absorption coefficient a(410) and a(670), as well seen in Figure 3. However, we can note that the ratio r rs (410)/r rs (670) is less under a CB than in its absence. We can explain that by the effect of spectral dependency of b b (λ) on the parameter u = b b /( a+ b b ); as seen from Figure 3, the b b coefficient is noticeably higher at 410 nm than at 670 nm for N coc = 10 7 cell/L. We can explain the observed changes by the behavior of the ratio of the backscattering and absorption coefficients b b /a that is seen in Figure 3.

Coccolithophore Bloom Impact on the Water Radiance Reflectance
In Section 3.1, we focused on the angular distributions of the underwater radiance including the distribution just beneath the sea surface that determines the so-called water-leaving radiance just above the sea surface. The main effect of the coccolithophore bloom on the upwelling radiance is the sharp increase in its absolute values, but the processing algorithms for satellite ocean color data are based on the analysis of the spectral values measured by satellite sensors. Keeping this in mind, we consider in this section the spectral values of the water radiance reflectance ρ(λ). Figure 9 shows the spectral values of ρ(λ) obtained by numerical calculation (see Section 2). We consider the changes of ρ(λ) depending on the coccolithophore concentration N coc for the Black Sea (the upper row in Figure 9) and the Barents Sea (the lower row). For each sea, we consider two cases of the spectral absorption (see Figure 3A

Coccolithophore Bloom Impact on the Water Radiance Reflectance
In Section 3.1, we focused on the angular distributions of the underwater radiance including the distribution just beneath the sea surface that determines the so-called water-leaving radiance just above the sea surface. The main effect of the coccolithophore bloom on the upwelling radiance is the sharp increase in its absolute values, but the processing algorithms for satellite ocean color data are based on the analysis of the spectral values measured by satellite sensors. Keeping this in mind, we consider in this section the spectral values of the water radiance reflectance ρ(λ). Figure 9 shows the spectral values of ρ(λ) obtained by numerical calculation (see Section 2). We consider the changes of ρ(λ) depending on the coccolithophore concentration Ncoc for the Black Sea (the upper row in Figure 9) and the Barents Sea (the lower row). For each sea, we consider two cases of the spectral absorption (see Figures   As seen for the Black Sea, an increase in the solar zenith angle leads to a slight change in ρ(λ). For the Barents Sea, water radiance reflectance for θ0 = 60° is greater than for overcast. The difference is at a maximum for Ncoc = 0 and a(440) = 0.18 m −1 , but even here it does not exceed 6%.
Comparing Figures 9A,C, it can be noted that all the curves in the blue part of the spectrum are much flatter for the Barents Sea than for the Black Sea. This can be easily explained: the spectral As seen for the Black Sea, an increase in the solar zenith angle leads to a slight change in ρ(λ). For the Barents Sea, water radiance reflectance for θ 0 = 60 • is greater than for overcast. The difference is at a maximum for N coc = 0 and a(440) = 0.18 m −1 , but even here it does not exceed 6%.
Comparing Figure 9A,C, it can be noted that all the curves in the blue part of the spectrum are much flatter for the Barents Sea than for the Black Sea. This can be easily explained: the spectral dependence of ρ is largely determined by the parameter u = b b /(a + b b ), and therefore by the coefficient a. Spectral dependences of ρ shown in the right panel of the figure are much more similar due to coincidence of absorption coefficients. For N coc = 10 7 cell/L, and to a large extent N coc = 5 × 10 6 cell/L, these curves are difficult to distinguish, since the scattering here is mainly determined by coccolithophores. At lower values of N coc , this similarity persists, although the scattering here is mainly due to non-coccolithophore suspension. The fact is that in this range of N coc values, the shapes of the ρ(λ) curves are largely determined by the spectral dependence of the absorption coefficient.
To track changes in the shape of the spectra ρ(λ) depending on the N coc concentration, we took the ratio ρ(490)/ρ(550) of the ρ spectral values at 490 and 550 nm; that is close to the parameters which are used in the processing algorithms OC2M and KD2M for the MODIS data [53]. Table 4 presents the results obtained. The data in Table 4 correspond to the changes in the spectra ρ(λ) in Figure 9. In the case of significant absorption (a(440) = 0.53 m −1 and 0.18 m −1 ), the ratio ρ(490)/ρ(550) in Table 4 depends weakly on N coc concentration and varies no more than 6%. For weak absorption (a(440) = 0.09 m −1 ), we observe a significant decrease in the ratio ρ(490)/ρ(550) with increasing N coc . Even for N coc = 10 6 cell/L, the value of ρ(490)/ρ(550) reduces by 10%, as compared with the bloom absence (N coc = 0); for an intense bloom (N coc ≥ 5 × 10 6 cell/L), the ratio ρ(490)/ρ(550) decreases by more than 26%. These results indicate the possibility of using the above parameter in the weak absorption case, and the need to choose another one for significant absorption.

Coccolithophore Impact on the Diffuse Attenuation Coefficient
In this section, we consider the propagation of the solar radiation in the water column, focusing on the effect of the coccolithophore bloom on this process. Usually, we assume for the attenuation of the underwater irradiance with depth an exponential law with an exponent K d called the diffuse attenuation coefficient. K d depends strongly on the wavelength and the angular structure of the propagating light stream; this problem was discussed in detail in [54]; We consider the applicability of the approximate formulas in Section 4.2. Figure 10 shows the vertical profiles of the downward irradiance E d (500) calculated for the case of the Black Sea, depending on the coccolithophore concentration N coc and the absorption coefficient a(440). As seen, both factors have a significant effect on the depth of the photic layer Z 1% (where PAR decreases to 1% of its surface value). In the case of a(440) = 0.09 m −1 and the CB absence, the Z 1% value is more than 40 m, it decreases to 34.5 m for N coc = 1 × 10 6 cell/L, and Z 1% is only 15.5 m for N coc = 10 7 cell/L. profiles of the downward irradiation on a semi-logarithmic scale ( Figure 10) are straight lines and there is no difference in the slope depending on the depth, but this is not so.   Figure 11 shows the difference between the profiles of the downward irradiation obtained as a result of numerical simulation, with the assumption that Kd is constant and equal to K0. In this case, we calculated the value K0 in the uppermost layer 0-0.1 m (as the increments of depth for calculations are 0.1 m) and then the irradiance profile using the exponential law with Kd = K0 = const. It is seen that even in the absence of a CB, the use of K0 overestimates the values of irradiance obtained; the difference for Ed estimates with a depth is especially large for small absorption and high coccolithophore concentrations.
Obviously, the use of K0 for calculating the attenuation of light with depth in the CB case results in large errors; it is better to use the Kd average for the photic layer. For example, in [54], the diffuse attenuation coefficient <K> was introduced as the average for layer τ10, in which the irradiance decreases by a factor of 10 compared to the surface value.  Table 5 presents the depth of the photic layer Z 1% for the cases of the Black and Barents Seas, depending on the coccolithophore concentration N coc , the absorption coefficient a(440), and illumination conditions. If the CB is absent, the depth of the photic layer Z 1% is much less in the Barents Sea than in the Black Sea, which is explained by a large amount of "non-coccolithophore" particles. In both seas, the intense CB leads to a significant reduction in the Z 1% value. For N coc = 10 7 cell/L, it decreases by a factor of two, regardless of the value of a(440). The Z 1% depends on the illumination conditions, in particular, on the low sun. The Z 1% is less than in overcast conditions, and in this case, it is less than with the high sun; the influence of illumination conditions mitigates with increasing N coc . Table 5. The depth of the photic layer Z 1% (m) for λ = 500 nm, in the cases of the Black Sea and Barents Sea, depending on the value of N coc , the illumination conditions, and the absorption coefficient.  It is known that the diffuse attenuation coefficient K d depends on the depth, even if the optical properties do not change with depth, but how can the CB affect these changes? It seems that the profiles of the downward irradiation on a semi-logarithmic scale ( Figure 10) are straight lines and there is no difference in the slope depending on the depth, but this is not so. Figure 11 shows the difference between the profiles of the downward irradiation obtained as a result of numerical simulation, with the assumption that K d is constant and equal to K 0 . In this case, we calculated the value K 0 in the uppermost layer 0-0.1 m (as the increments of depth for calculations are 0.1 m) and then the irradiance profile using the exponential law with K d = K 0 = const. It is seen that even in the absence of a CB, the use of K 0 overestimates the values of irradiance obtained; the difference for E d estimates with a depth is especially large for small absorption and high coccolithophore concentrations. The values of the spectral coefficient of diffuse attenuation <K> obtained as a result of the simulation are shown on Figure 12. For the Black Sea, calculations were carried out for θ0 = 25°, and for the Barents Sea, -θ0 = 60°. It is seen that an increase in Ncoc leads to a significant increase in <K>, especially in the case of low absorption. For example, for the Black Sea, λ = 500 nm, and a(440) = 0.09 m -1 , the diffuse attenuation coefficient <K> equals 0.094 m −1 in the absence of a CB, it increases by a factor of two for Ncoc = 3 10 6 cell/L (0.183 m −1 ) and by three for Ncoc = 10 7 cell/L (0.304 m −1 ); for the Barents Sea and a(440) = 0.18 m −1 , the diffuse attenuation coefficient <K> equals 0.28 m −1 in the absence of a CB and it increases only by a factor of two for Ncoc = 10 7 cells/l (0.56 m −1 ), unlike the case of Black Sea. For a(440) = 0.53 m −1 , the increase in <K> does not occur so rapidly with increasing Ncoc. Due to a large amount of "non-coccolithophore" particles in the Barents Sea, the obtained <K> values are significantly higher than similar estimates for the Black Sea (see Figure 12 and Table 6). Table 6 shows the results of the <K> dependence from the illumination conditions for the Black Sea and Barents Sea. For both seas, the highest value of <K> is obtained in the case of a clear sky with a low sun (θ0 = 60°) and the smallest one with a high sun (θ0 = 25°).  Obviously, the use of K 0 for calculating the attenuation of light with depth in the CB case results in large errors; it is better to use the K d average for the photic layer. For example, in [54], the diffuse attenuation coefficient <K> was introduced as the average for layer τ 10 , in which the irradiance decreases by a factor of 10 compared to the surface value.
The values of the spectral coefficient of diffuse attenuation <K> obtained as a result of the simulation are shown on Figure 12. For the Black Sea, calculations were carried out for θ 0 = 25 • , and for the Barents Sea, -θ 0 = 60 • . It is seen that an increase in N coc leads to a significant increase in <K>, especially in the case of low absorption. For example, for the Black Sea, λ = 500 nm, and a(440) = 0.09 m −1 , the diffuse attenuation coefficient <K> equals 0.094 m −1 in the absence of a CB, it increases by a factor of two for N coc = 3 × 10 6 cell/L (0.183 m −1 ) and by three for N coc = 10 7 cell/L (0.304 m −1 ); for the Barents Sea and a(440) = 0.18 m −1 , the diffuse attenuation coefficient <K> equals 0.28 m −1 in the absence of a CB and it increases only by a factor of two for N coc = 10 7 cells/l (0.56 m −1 ), unlike the case of Black Sea. For a(440) = 0.53 m −1 , the increase in <K> does not occur so rapidly with increasing N coc . Due to a large amount of "non-coccolithophore" particles in the Barents Sea, the obtained <K> values are significantly higher than similar estimates for the Black Sea (see Figure 12 and Table 6). Table 6 shows the results of the <K> dependence from the illumination conditions for the Black Sea and Barents Sea. For both seas, the highest value of <K> is obtained in the case of a clear sky with a low sun (θ 0 = 60 • ) and the smallest one with a high sun (θ 0 = 25 • ).
Lee et al. [59] modified (9) for both the coastal and open ocean waters, suggesting to use it with other coefficients ρ(u) = π (0.089 + 0.1245 u) u (10) Figure 13 shows the plots of ρ(λ), according to Equations (7)-(10) (a similar figure for Equations (7)-(9) was presented in [56]). We see that all Equations (7)-(10) are very close to each other for values u ≤ 0.25 when ρ(λ) ≤ 0.1, a noticeable difference appears only for u ≥ 0.2. For large values of the parameter u, corresponding to the high coccolithophore concentration, Equation (7) gives incredibly large values of ρ; Equations (8) and (10) give similar results; Equation (9) gives the lowest value of ρ, as compared to other formulas. (10) Figure 13 shows the plots of ρ(λ), according to Equations (7-10) (a similar figure for Equations (7-9) was presented in [56]). We see that all Equations (7-10) are very close to each other for values u ≤ 0.25 when ρ(λ) ≤ 0.1, a noticeable difference appears only for u ≥ 0.2. For large values of the parameter u, corresponding to the high coccolithophore concentration, Equation (7) gives incredibly large values of ρ; Equations (8) and (10) give similar results; Equation (9) gives the lowest value of ρ, as compared to other formulas. To assess the applicability of Equations (7-10) under the coccolithophore bloom conditions, a numerical simulation of the dependence of ρ on the coccolithophore concentration was made. We have set the acceptable absolute error for ρ equal to 0.01 for ρ ≤ 0.1 (u ≤ 0.25), and the acceptable relative error equal to 10% for ρ > 0.1 or u > 0.25.
It should be noted that all Equations (7-10) do not depend on the solar zenith angle θ0. As shown in Section 3.2, a difference between the values of ρ in various atmospheric conditions is quite insignificant. Keeping it in mind, we performed our numerical calculations, intended for the optimal choice between Equations (7-10) for the cloudy sky with an optical thickness of cloud τcl = 20, since it gives the value of ρ close to the average for the three considered cases of illumination conditions To assess the applicability of Equations (7)-(10) under the coccolithophore bloom conditions, a numerical simulation of the dependence of ρ on the coccolithophore concentration was made. We have set the acceptable absolute error for ρ equal to 0.01 for ρ ≤ 0.1 (u ≤ 0.25), and the acceptable relative error equal to 10% for ρ > 0.1 or u > 0.25.
It should be noted that all Equations (7)-(10) do not depend on the solar zenith angle θ 0 . As shown in Section 3.2, a difference between the values of ρ in various atmospheric conditions is quite insignificant. Keeping it in mind, we performed our numerical calculations, intended for the optimal choice between Equations (7)-(10) for the cloudy sky with an optical thickness of cloud τ cl = 20, since it gives the value of ρ close to the average for the three considered cases of illumination conditions (clear sky with θ 0 = 25 • and 60 • , and overcast). Figure 14 presents the relative errors of Equations (7)-(10) as a function of the parameter u for u ≥ 0.25, various values of the absorption coefficient a, taking into account "non-coccolithophore" scattering particles, and without it. We estimated the relative errors of Equations (7)-(10) at 530 nm, performing our numerical simulation for N coc from 0 to 10 7 cell/L. The two upper pictures in Figure 14  Relative errors of Equations (7-10), depending on the parameter u, taking into account the "non-coccolithophore" particles (dashed lines) and without it (solid lines). The color of the curves is the same as in Figure 13.
It is seen that Equation (8) and Equation (10) result in the acceptable error for all cases under consideration, in contrast to Equation (7) and Equation (9). We explain it by the fact that the Equation (8) and Equation (10) were derived for highly scattering waters. For u < 0.5, the Equation (8) is usually better than (10), for u > 0.5 Equation (10) it is otherwise. It is well seen for the highest u values in the Black Sea in the case of low absorption ( Figure 14C). In the case of u < 0.25, all formulas result in 0.01, and all of them are acceptable. Figure 15 shows the results of comparisons between the spectral values of ρ(λ) obtained by numerical simulation and the estimates by Equations (7-10) for the Barents Sea and the Black Sea. In the Black Sea, the numerical simulations were performed with a clear sky and θ0 = 25°, while in the Barents Sea at θ0 = 60°. We calculated the ρ(λ) values for Ncoc = 10 7 cell/L and low absorption, corresponding to the maximum u values. It is seen in Figure 15 that the values of ρ(λ) by the numerical calculations and Equations (7-10) are very close to each other in the case of u ≤ 0.25 for the wavelength λ > 600 nm, the lowest approximation errors are obtained by Equations (8) and (10) for u > 0.25 and Figure 14. Relative errors of Equations (7)-(10), depending on the parameter u, taking into account the "non-coccolithophore" particles (dashed lines) and without it (solid lines). The color of the curves is the same as in Figure 13.
The values of relative error are shown at the vertical axis (ρ appr −ρ)/ρ, where ρ is the result of the numeric simulation and ρ appr is calculated by one of the Equations (7)-(10). We took into account the "non-coccolithophore" particles (shown by dashed lines) and the solid lines show the relative errors without "non-coccolithophore" particles. In Figure 14, the minimum value of u = 0.25 corresponds ρ = 0.1, as we set the acceptable relative error (10%) for ρ > 0.1. The maximum values of the parameter u were obtained for N coc = 10 7 cell/L, and they are much larger in the case of low absorption.
It is seen that Equations (8) and (10) result in the acceptable error for all cases under consideration, in contrast to Equations (7) and (9). We explain it by the fact that the Equations (8) and (10) were derived for highly scattering waters. For u < 0.5, the Equation (8) is usually better than (10), for u > 0.5 Equation (10) it is otherwise. It is well seen for the highest u values in the Black Sea in the case of low absorption ( Figure 14C). In the case of u < 0.25, all formulas result in 0.01, and all of them are acceptable. Figure 15 shows the results of comparisons between the spectral values of ρ(λ) obtained by numerical simulation and the estimates by Equations (7)-(10) for the Barents Sea and the Black Sea. In the Black Sea, the numerical simulations were performed with a clear sky and θ 0 = 25 • , while in the Barents Sea at θ 0 = 60 • . We calculated the ρ(λ) values for N coc = 10 7 cell/L and low absorption, corresponding to the maximum u values. It is seen in Figure 15 that the values of ρ(λ) by the numerical calculations and Equations (7)-(10) are very close to each other in the case of u ≤ 0.25 for the wavelength λ > 600 nm, the lowest approximation errors are obtained by Equations (8) and (10) for u > 0.25 and λ < 600 nm, and the best approximation for u > 0.5 in the case of the Black Sea, is obtained by Equation (10).   In the presence of a CB, even one that is not so intensive (Ncoc = 1 × 10 6 cell/L, u = 0.34), the influence of the solar zenith angle weakens and acceptable values of the relative errors are obtained for the Equations (7-8, 10). For Ncoc ≥ 3 and 10 6 cell/L, u > 0.5 a good approximation and is obtained only by the Equation (8) and Equation (10), regardless of the solar zenith angle.  Table 7 presents the absolute and relative errors of Equations (7)-(10) depending on the value of N coc and the solar zenith angle for the Black Sea, a(440) = 0.09 m −1 and λ = 490 nm, where the parameter u reaches its maximum value. In the absence of a CB for u = 0.17, the errors of the formulas increase significantly for the large solar zenith angle, but for all formulas, the absolute error does not exceed the threshold value of 0.01.
In the presence of a CB, even one that is not so intensive (N coc = 1 × 10 6 cell/L, u = 0.34), the influence of the solar zenith angle weakens and acceptable values of the relative errors are obtained for the Equations (7), (8) and (10). For N coc ≥ 3 and 10 6 cell/L, u > 0.5 a good approximation and is obtained only by the Equations (8) and (10), regardless of the solar zenith angle.
Summing up the results of our calculations, we can conclude that for u ≤ 0.25, any Equation of (7)-(10) allows us to estimate ρ with an accuracy not worse than 0.01. For u > 0.25, the Equations (8) or (10) give an accuracy not worse than 10%; for u > 0.5, Equation (10) is more reliable.

Diffuse Attenuation Coefficient
Gordon [54] considered two variants for the diffuse attenuation coefficient K d . The first case is K 0 just beneath the sea surface, the second, <K>, averaged over half of the euphotic zone, corresponding to the optical depth τ 10 = 2.3 [54]. In Section 3.3, we calculated both coefficients by using our numerical model and estimated the errors in the downward irradiance E d (z) obtained by using the exponential law with the constant diffuse attenuation coefficient. In this section, we are estimating the possibility of the use of the approximate formulas for K d , two of them are Gordon's formulas [54], the third is developed by Lee et al. [60].
First, Gordon's formula for K 0 : where D o is the parameter of the angular radiance distribution just beneath the surface for a fully absorbing ocean. D o depends on the solar zenith angle θ 0 , atmospheric conditions (clear or overcast sky), and the wavelength (some results of the computed D o are given in Table 3 [54]). This simple formula is the most often used to calculate K d through the inherent optical characteristics of seawater. Gordon's formula for <K> is a polynomial of the third degree from (1−ω 0 F): < k > n a + bb c n (12) where F is the forward scattering probability, ω 0 is the single scattering albedo, and c is the beam attenuation coefficient; the <k> n coefficients were derived from the results of simulations and its values are given in [54]. For applying Equation (12) The formula was derived from the simulations with a clear sky, so it cannot be used for cloudy conditions. Figure 16 and Table 8 present the accuracy estimates for the above formulas depending on the CB intensity in the Black Sea. Figure 16A,B show the applicability of Equation (11) for the estimation of K 0 . As seen in this case, the errors are practically independent of the CB presence; the average error equals 0.01-0.04 m −1 , the maximum is 0.1 m −1 . Figure 16C,D demonstrate that the errors of Equation (12) are noticeably larger for N coc = 10 7 cell/L and weakly depend on the absorption spectrum. For an intense CB, Equation (12)    As mentioned above, for estimating the attenuation of underwater irradiance with depth, a simpler Equation (11) is often used instead of <K>. Figure 16E,F show that it is acceptable in the absence of a CB, but this formula results in noticeably increasing errors in the case of an intense CB. In this case, Equation (11) underestimates K d by 0.06-0.21 m −1 (18-33%) for a small a and by 0.17-0.43 m −1 (20-33%) for a large a. Figure 16G,H show the results of calculations by Equation (13), as compared with numerical calculations. In the absence of a CB, the error of this formula does not exceed 0.03 m −1 (4%), but in the case of N coc = 10 7 cell/L and the absorption value of a(440) = 0.53 m −1 , Equation (11) Table 8 presents the absolute errors of the Equations (11)-(13) depending on the N coc value, the solar zenith angle θ 0 , and the absorption coefficient a(440). As seen, all formulas give good results in the absence of a CB. For higher N coc values, the errors increase, and only Equations (12) and (13) formulas can provide reasonable results. We can see that in most cases, Equation (12) gives better results for a(440) = 0.09, (11) for a(440) = 0.53. For θ 0 = 60 • , the Equation (13), including the solar zenith angle θ 0 as the input parameter, gives acceptable results for within the whole spectral range for N coc = 5·10 6 cell/L and for N coc = 10 7 only for a(440) = 0.53 (for a(440) = 0.09 only for the range of 600-700 nm).
The test results of the accuracy of Equations (11)-(13) for the Barents Sea are similar to the ones obtained for the Black Sea. Figure 17 presents, as an example, the results of comparing the <K> spectra in the Barents Sea (θ 0 = 60 • ) obtained by numerical simulation and by Equations (12) and (13), depending on the value of N coc and the absorption coefficient a(440). As in the Black Sea, Equation (12) underestimates the value of <K> for an intense CB, but by a smaller value. For N coc = 10 7 cell/L Equation (12) underestimates <K> by no more than 0.08 m −1 when using the absorption spectrum with a(440) = 0.18 m −1 measured by ICAM and by no more than 0.14 m −1 in the case of the model absorption spectrum with a(440) = 0.53 m −1 . As in Table 8, the errors of the Equation (13)  Thus, in the CB case, a choice between Equations (12) and (13) depends on the absorption coefficient value. If the bb value is not too high, compared with the absorption value, it is better to use a simple Equation (13). But this formula results in large errors in the case of low absorption and intense scattering (bb ≥ 0.5 a); in this case, it is better to use Equation (12). However, for its application, knowledge of the bb/b ratio is needed. This raises a problem in its application in the case of remote sensing. The use of Equation (11) is possible only in the CB absence.
For the estimation of the diffuse attenuation coefficient value in the absence of a CB under clear sky conditions, you can use any of Equations (11)(12)(13); in the case of the cloudy sky, Equations (11) or (12). In the CB case and a clear sky, if the absorption is not too low (bb ≤ 0.5 a), you can use Equation (13). In the case of high scattering and low absorption, you can use Equation (12), if you know the value of bb/b. Thus, in the CB case, a choice between Equations (12) and (13) depends on the absorption coefficient value. If the b b value is not too high, compared with the absorption value, it is better to use a simple Equation (13). But this formula results in large errors in the case of low absorption and intense scattering (b b ≥ 0.5 a); in this case, it is better to use Equation (12). However, for its application, knowledge of the b b /b ratio is needed. This raises a problem in its application in the case of remote sensing. The use of Equation (11) is possible only in the CB absence.
For the estimation of the diffuse attenuation coefficient value in the absence of a CB under clear sky conditions, you can use any of Equations (11)- (13); in the case of the cloudy sky, Equations (11) or (12). In the CB case and a clear sky, if the absorption is not too low (b b ≤ 0.5 a), you can use Equation (13). In the case of high scattering and low absorption, you can use Equation (12), if you know the value of b b /b.

Discussion and Conclusions
In this section, we will focus on the key issues of our modeling, the choice of a calculation method, forming a set of input parameters, a discussion of one of the most interesting results-the asymptotic regime-and formulate the main conclusions, in particular, the scientific and practical significance of the results.

Comments to the Calculation Methods
In the Introduction, we considered the most widely used numerical methods: Monte Carlo, DISORT, and Hydrolight. The last seems to be the most appropriate, having a high processing speed and accuracy comparable to Monte Carlo (M-C). But this is a commercial program, rather expensive, and we do not have it. The Monte Carlo method works slowly, but a good agreement with Monte Carlo indicates the reliability of the results. Therefore, we used the modified DISORT for our numerical calculations (see Section 2.2) and M-C for checking the obtained results in selected cases. As an example, we present the results of such a comparison between the underwater radiance calculation by Monte Carlo-and DISORT-based matrix operator methods ( Figure 18).

Discussion and Conclusion
In this section, we will focus on the key issues of our modeling, the choice of a calculation method, forming a set of input parameters, a discussion of one of the most interesting results-the asymptotic regime-and formulate the main conclusions, in particular, the scientific and practical significance of the results.

Comments to the Calculation Methods
In the Introduction, we considered the most widely used numerical methods: Monte Carlo, DISORT, and Hydrolight. The last seems to be the most appropriate, having a high processing speed and accuracy comparable to Monte Carlo (M-C). But this is a commercial program, rather expensive, and we do not have it. The Monte Carlo method works slowly, but a good agreement with Monte Carlo indicates the reliability of the results. Therefore, we used the modified DISORT for our numerical calculations (see Section 2.2) and M-C for checking the obtained results in selected cases.
As an example, we present the results of such a comparison between the underwater radiance calculation by Monte Carlo-and DISORT-based matrix operator methods ( Figure 18).

Set of the Input Parameters
We indicated in Section 2.1 that we did not set out for our modeling to reproduce the results of the field measurements. Still, we needed in situ measured data to approximate the input modeling parameters to the real conditions observed in the Barents Sea and Black Sea. We did not measure all input parameters needed for the calculations, but the measured data allowed us to derive the unmeasured ones (see Table 1).
The choice of appropriate optical characteristics for the main CB components-plated cells and detached coccolites-is particularly difficult due to their complicated shape and noticeable variability depending on external conditions, in particular, on temperature and mineral nutrition. Balch [52] noted "changes in coccolith morphology would cause variability in their backscattering crosssections, which would translate to changes in their reflectance." Von Dassow et al. [61] reported about the coccoliths' ability to rotate linearly polarized light.
Interest in optical characteristics has increased in the last twenty years in connection with the

Set of the Input Parameters
We indicated in Section 2.1 that we did not set out for our modeling to reproduce the results of the field measurements. Still, we needed in situ measured data to approximate the input modeling parameters to the real conditions observed in the Barents Sea and Black Sea. We did not measure all input parameters needed for the calculations, but the measured data allowed us to derive the unmeasured ones (see Table 1).
The choice of appropriate optical characteristics for the main CB components-plated cells and detached coccolites-is particularly difficult due to their complicated shape and noticeable variability depending on external conditions, in particular, on temperature and mineral nutrition. Balch [52] noted "changes in coccolith morphology would cause variability in their backscattering cross-sections, which would translate to changes in their reflectance." Von Dassow et al. [61] reported about the coccoliths' ability to rotate linearly polarized light.
Interest in optical characteristics has increased in the last twenty years in connection with the task of quantifying CB parameters from satellite color scanners. In addition to the works [47][48][49][50][51][52] listed in Section 2.3, there are several other publications in which optical parameters were estimated theoretically or from field measurements [62][63][64]. However, the representative set of optical parameters of CB components is only given in the publication [51]. That is the main reason why we took [51] as the base for our numerical calculations. The other reason is the possibility to use the data of [51] for the quantitative assessment of the intensity of coccolithophore blooms in the Barents Sea from satellite data [44].

Asymptotic Regime
Many researchers have studied the asymptotic regime and derived formulas to describe its angular radiance distribution and calculate the diffuse attenuation coefficient; the interested reader can find the references, including works from the first decade of this millennium, in the books of Kirk [13] and Gordon [21]. All authors agree that the single-scattering albedo, the single-scattering average cosine, and the average square angle are the main parameters determining the features of the asymptotic regime.
We pay attention to the theoretical results presented in [18]. The approach, based on solving the characteristic equation in the small-angle diffusion approximation, allowed the derivation of the formulas for the asymptotic radiance distribution in the range of angles 0 • -45 • and for estimating the optical depth corresponding to the establishment of the asymptotic regime [18,65]. The examples of estimating such depths for the clearest (1) and turbid coastal (2) waters are presented in [18].
The following parameters were taken at 555 nm: (1) The beam attenuation coefficient c(555) = 0.07 m −1 , the single scattering albedo ω 0 (555) = 0.4, the scattering coefficient b(555) = 0.028 m −1 , the average single-scattering cosine <cosθ> = 0.8; (2) 1.0 m −1 , 0.9, 0.9 m −1 , and 0.48, respectively. The boundary estimates are (1) 23 and 330 m, (2) 60 and 60 m, for the optical and geometrical depth, respectively. It would be interesting to compare these estimates with those obtained from our calculations. As mentioned above, two main conditions determine the asymptotic regime: the angular radiance distribution is symmetrical concerning the vertical axis and unchangeable with depth; The diffuse attenuation coefficient is constant over depth. Table 3 shows the necessary parameters for the depth estimation.
Because we deal with the asymptotic regime, we should set requirements for the accuracy of the fulfillment of the above conditions. As we noted in the Introduction, the advantage of numerical calculations is that their accuracy can be much higher than practically possible in field measurements. We set an accuracy of 3% for the angular distributions and 0.01 m −1 for the diffuse attenuation coefficient k. As can be seen from the Table 3, taking into account the depth discreteness of the presented data, we can accept the following values of the optical and geometrical depth as the boundary of the onset of the asymptotic regime: 6 m and 8 m, respectively, for N coc = 5 × 10 6 cell/L, a(440) = 0.09 and 0.53 m −1 , and 5 m for both cases of N coc = 10 7 cell/L. The values of the optical depth for the above boundaries equal 17, 24, 26, and 26, respectively. As compared with the results given in [18], we can note the close values of the optical depths for their clearest case, but the geometrical depths differ by factors of 40-50. This is an important result from several points of view which have practical significance. In reality, we have a light curtain that creates a significant obstacle for observing underwater objects at shallow depths, as well as for transmitting signals across the surface when using laser communication with underwater objects. On the other hand, we get a natural laboratory, allowing us to study the features of the asymptotic regime in the subsurface layer; the cost of such research, no doubt, is much less than deep-sea research. Such studies can provide a possibility to develop fairly simple methods for evaluating a set of inherent optical characteristics, including the parameters of the phase scattering function.

Conclusions
The coccolithophore blooms are a planetary-scale phenomenon, a potent producer of CaCO 3 in the ocean, affecting the exchange of carbon dioxide between the ocean and atmosphere, the greenhouse effect, and global climate change significantly. We studied the impact of intense coccolithophore blooms on the underwater light fields, which provide the most comprehensive information about the propagation of solar radiation in the water column. We prepared software and the sets of input parameters, corresponding to real conditions in the Barents Sea and Black Sea, for modeling the underwater light fields and their basic parameters under extreme coccolithophore bloom conditions. We used the modified DISORT code and the Monte Carlo method to get the exact data of our numerical calculations. The most interesting of the obtained results comprise data on the asymptotic regime features, whose border of onset can be about five meters as compared with a few hundred meters in the open ocean. This result is of practical significance, denoting an obstacle for observing underwater objects at shallow depths and transmitting the optical signals across the sea surface. On the other hand, we get a natural laboratory, allowing us to study the features of the asymptotic regime near the sea surface. The assessment of the accuracy of the approximating formulas with our exact data enables us to make recommendations for the choice of the optimal algorithm under intensive coccolithophore bloom conditions.