1. Introduction

Scattering is a fundamental optical process which largely governs radiative transport processes within the visible spectrum in the ocean and the atmosphere [1]. In optical oceanography, measuring scattering properties of the water column can give valuable information about size and composition of particles submerged in the water mass, such as phytoplankton or colloidal mineral particles [2–9]. Both the scattering and backscattering coefficient have been found to have a strong relationship with suspended particulate matter concentration [2–4]. These properties tell us how much light is scattered in the backward direction (backscattering) and in total (scattering). The ratio between the backscattering and scattering coefficient yields an estimate of the bulk refractive index, which indicates whether inorganic or organic particles are dominating [5]. The full angular distribution of scattering, namely the volume scattering function (VSF), has not been widely measured due to difficulties with designing instrumentation and conducting measurements in situ, but may tell us much more about the particulate properties [6,7]. The VSF at angles smaller than 15$^{\circ }$ is now routinely used for estimating the particle size distribution between approximately 1 and 500 $\mu$m [8,9]. The VSF is also an important property to parameterize correctly in radiative transfer modelling of the ocean, which is critical for achieving accurate ocean color remote sensing observations of the oceanic and coastal environments around the globe [10,11]. Ecologically and economically important coastal environments are of considerable interest to measure with in situ instrumentation, as these environments are still difficult to observe accurately using satellite observations due to optical and spatio-temporal complexity [12]. Scattering properties in coastal and estuarine environments have a large dynamic range, for instance, the scattering coefficient can span from approximately 0.05 to 30 m$^{-1}$ [2], which puts extra demands on instrumentation.

The LISST-VSF (Sequoia Scientific) is a commercially developed instrument able to measure the volume scattering function (VSF) from ${0.1}^{\circ}$ to ${150}^{\circ}$. From the VSF one can also derive other important inherent optical properties such as the scattering coefficient and phase function. Recent studies [13,14] have demonstrated that, for sufficiently large optical depths, multiple scattering results in significant errors in the VSF measured by the LISST-VSF. This error arises from the assumption that all detected light is only scattered once. The detected signal is corrected for light lost, either to scattering or absorption, along the path from laser to detector by multiplying the detected intensity with the factor $e^{cl}$. Here, $c$ is the attenuation coefficient and $l$ is the path length from laser to detector. However, not all detected light is scattered only once. Some of the light which is assumed to be lost due to scattering along the path, may still reach a detector. Thus, these photons add to the single-scattering signal, resulting in elevated VSFs, scattering coefficients $b$, and backscattering coefficients $b_b$, as well as distorted phase functions.

Monte Carlo simulations have been widely used to investigate measurement errors in optical instruments, as it is relatively easy to implement complex geometries without significant approximations [15–21]. For instance, scattering errors in spectrophotometric absorption meters have been assessed using Monte Carlo methods in several studies including Kirk [15], Stramski and Piskozub [16], and McKee et al. [17]. Piskozub et al. [18] found that multiple scattering has a negligible effect on attenuation measurements. Some recent research efforts with Monte Carlo modelling have been directed towards backscattering sensors. Doxaran et al. [19] simulated two backscattering sensors and found that the ECO-BB instrument (Sea-Bird Sci) is less affected by scattering errors than HydroScat instruments (HOBI Labs), due to a shorter path length (0.02 m, compared to 0.10 m). Vadakke-Chanat et al. [20] also conducted an assessment of the ECO-BB instrument with Monte Carlo modelling and reproduced many of the findings in Doxaran et al. [19], but the effective path length found for the instrument differs between both studies and the user manual. Recently, another study combining direct modeling of light propagation and Monte Carlo simulations of the HydroScat and ECO-BB instruments was published [21]. Here, angular response functions were derived for both instruments, and ECO-BB was found to have an especially large angular full width at half maximum, which can be a significant error source, while multiple scattering due to long path lengths is the main uncertainty for the HydroScat instruments. Monte Carlo modelling has also been applied to other aspects of marine and atmospheric optics, including radiative transfer, underwater optical wireless communication, and LIDAR [22–24].

The effect of multiple scattering on VSF instruments has been studied less than the error for backscattering sensors. In a recent study [14], we presented a Monte Carlo simulation developed to analyze and quantify multiple scattering errors in the LISST-VSF. The simulation was compared to experimental measurements using ${508}\;\textrm{nm}$ polystyrene beads and Arizona test dust as the scattering agents. The simulation was found to accurately predict experimental measurements for the samples investigated. When comparing simulation to experimental results, the relative error did not exceed 13% for the scattering coefficient and 10% for the phase function for optical depths < 4.9, where a considerable fraction of these errors can be attributed to errors in the experimental measurements, e.g. detector saturation. Comparing simulated and theoretical VSF, i.e the relative measurement error originating from multiple scattering, the error was found to depend on the phase function of the sample volume. The relative error in the scattering coefficient reached 100% at an optical depth of 2.8 and 2.0 for the 508 nm polystyrene beads and Arizona test dust samples, respectively.

While the simulation is proven accurate, it suffers from long computation time due to the fact that most phase functions in natural waters are very forward peaked. Thus, a significant number of photons must be simulated in order to achieve low variance results in the backward direction. In this study, we present three variance reduction methods (VRMs) that have been implemented to increase the efficiency of the simulation. The VRMs are called mean free path transform (MFPTR), detector directional importance sampling (DDIS) and splitting, and are presented in section 2.2. Additionally, the simulation is used to further investigate the effects of multiple scattering on LISST-VSF measurements, especially for realistic phase functions. By using a range of different Fournier-Forand (FF) phase functions [25,26], we analyze how the errors originating from multiple scattering depend on the phase function, scattering coefficient, and absorption coefficient of the sample volume. This is presented in section 3.2.

2. Methods

2.1 Monte Carlo simulation

The Monte Carlo simulation presented in [14] is a straight forward algorithm where the individual photons are tracked through the sample volume until they cross the sample volume border or reach a detector. There are two separate detectors in the LISST-VSF instrument. The geometry of the simulation is presented in Fig. 1, and a more detailed description of the simulation geometry is given in [14]. A set of ring detectors are located at the bottom of the sample volume chamber, detecting light at scattering angles 0.01-${15}^{\circ}$, while scattering angles 15-${150}^{\circ}$ are detected by a rotating eyeball detector. The ring detectors are simulated as a circular surface, where photons crossing this surface are detected and sorted according to their angle of incidence. The eyeball detector is simulated as a ring torus, where crossing the surface area corresponds to detection by the eyeball. The path length between each scattering event is sampled from a distribution derived from the probability density function (PDF) calculated from the beam attenuation $f(l)$, and the scattering angle $\theta$ is sampled from a distribution derived from a PDF calculated from the phase function $f(\theta )$.

 

Fig. 1. Simulation geometry. The sample volume is limited to the cylinder. The detector areas are colored in blue, where the disk represents the set of ring detectors, and the ring torus represents the eyeball detector. The solid green arrow represent incident laser beam, while the dashed green arrows represent scattered light.

Download Full Size | PDF

2.2 Variance reduction methods

The VRMs presented here are based on manipulating the natural scattering behaviour of the photons. Photons are given a statistical weight, initially $w_0=1$, which can be adjusted in order to reflect the alteration in the natural scattering behavior. A common example of this is to account for absorption of the sample volume by scaling the weight at each scattering event to reflect the probability of absorption. In this case, the weight after a scattering event is calculated as $w_{i+1}=w_i \cdot b/c$, where $w_i$ is the weight before the scattering event. The first VRM presented in this study is called mean free path transform (MFPTR). Usually, this method is used to stretch the path length in order to probe deeper into a medium (e.g. in Mainegra-Hing and Kawrakow [27]). Here, the method is applied to limit the path length such that all photons are scattered within the sample volume, i.e. the transmission is zero. Thus, all photons interact with the sample volume, increasing the simulated signal at the detectors. The second VRM presented, addresses the issue of poor statistics in the backwards direction. The method is based on the work presented in a study by Buras and Mayer [22], and is called detector directional importance sampling (DDIS). In essence, DDIS increases the number of photons scattered towards the detector by giving each photon a probability of being turned directly towards the detector prior to sampling a scattering angle. A consequence of this method is the possibility of photons accumulating weight, which can result in spikes in the detected signal. Thus, the third method presented in this study is implemented to avoid large photon weights and is called splitting [22,27,28]. When the photon weight becomes larger than a threshold value, the photon is split into $N$ identical parts which are propagated individually from the position of splitting. Thus, the weight is distributed to more than one detection angle, avoiding spikes in the signal.

2.2.1 Mean free path transform (MFPTR)

Depending on the scattering coefficient of the sample in question, the fraction of transmitted photons can be large, meaning that a lot of simulation time is spent initiating photons that do not contribute to the detected signal. Here, we have applied MFPTR to the path length distribution before the first scattering event so that all photons are scattered within the sample volume. This is achieved by defining the path length PDF as

2.2.2 Detector directional importance sampling (DDIS)

In general, the signal is much higher for the ring detector than for the eyeball detector. Thus, DDIS is only applied to the eyeball detector, i.e. the probability of being scattered towards the eyeball detector is increased and the probability of being scattered towards the ring detector is decreased. Due to the complexity of the detector geometry in this simulation, the implementation of DDIS is not straight forward. First, a random azimuth angle $\phi$ is sampled from a uniform distribution [0, 2$\pi$]. Once an azimuth angle has been sampled, all possible scattering directions lies in a plane defined by the sampled $\phi$. By determining where this plane intersects the eyeball detector (ring torus), we find the angles at which the photon is scattered directly towards the eyeball detector. This is illustrated in Fig. 2(a) for a photon at scattering position P and direction prior to scattering indicated by the red arrow. The angles between the red arrow and the two yellow arrows in Fig. 2(a) indicate the scattering angles $\theta _{\textrm {rot},1}$ and $\theta _{\textrm {rot},2}$ giving two possible scattering directions directly towards the detector. Now, instead of letting it follow the natural path, the photon is given a probability $\epsilon$ of being redirected towards the eyeball detector through either the angle $\theta _{\textrm {rot},1}$ or $\theta _{\textrm {rot},2}$ (chosen with equal probability). Then the scattering angle $\theta _{\textrm {ddis}}$ is sampled according to a PDF $f_{\textrm {ddis}}(\theta _{\textrm {ddis}})$ tailored for the specific geometry in this simulation (see discussion below). Thus, the final scattering angle becomes

By contrast, with a probability of (1-$\epsilon$), the photon is scattered naturally, i.e. according to the phase function of the sample volume. In order to account for this alteration of the natural scattering behaviour, the weight of the photon must be scaled according to the function

Next, we want to choose $\Delta \theta$ so that we maximize the number of photons being scattered towards the detector, without compromising the total intensity (weight) detected. Since we know the distance $D$ from the photon position P to the eyeball detector, we can calculate the maximum angle at which the photon can scatter away from the center of the detector and still be within the detector area, see Fig. 2(b),

 

Fig. 2. Illustration of the geometry applied to find (a) the scattering direction resulting in the photon being sent directly towards the eyeball detector (blue ring torus), and (b) the maximum angle the photon can scatter away from the center of the detector while still being detected.

Download Full Size | PDF

2.2.3 Splitting

As mentioned, applying DDIS may result in photons accumulating weight, causing spikes in the detected signal. In order to avoid this, a photon is split into $N$ identical parts when the photon reaches a weight larger than a threshold value. When the photon is split, the weight of each part can be calculated as

The $N$ parts are then propagated individually from the position where the original photon was split and in the same scattering direction. For the results presented in Fig. 3, the threshold weight was chosen to be $w_{\textrm {max}}=1$, and $N$ was chosen to be 2 times the floor of $w$ ($N=2\lfloor {w}\rfloor$), so that the photon is split into a minimum of two parts. The threshold value can be set to $w_{\textrm {max}}=1$ because the MFPTR results in reduced weight before the first scattering event. However, this assumption may fail when the optical depth becomes very large.

2.3 Analysing simulated VSFs

Using the improved MC simulation described in section 2.2, a range of different sample volumes were simulated using FF phase functions with scattering coefficients $b$ ranging from 0.05 to ${20}\;\textrm{m}^{-1}$. The wide range of different phase functions used enables the investigation of the multiple-scattering errors dependency on the phase function. The FF phase function is an approximate analytic expression of the phase function for an ensemble of spherical particles [25,26]. It has been extensively used for modeling the VSF of natural waters, usually providing a better fit than the simple one-term Henyey-Greenstein phase function. The function takes two parameters $n$ and $m$, where $n$ is the refractive index of the particles and $m$ defines the particle size distribution according to the inverse power-law $f(D)\propto D^{-m}$, where $D$ is particle diameter. The quantitative investigation of the multiple-scattering errors is performed by calculating the relative error for the scattering coefficient $b$, backscattering coefficient $b_b$ and phase function $p(\theta )$. The relative error in the scattering coefficient is calculated as

Thus, $b_{\textrm {true}}$ is in fact slightly smaller than the true scattering coefficient, but equal to the scattering coefficient one would measure given a perfect measurement, i.e. no multiple scattering. The relative error $e_{b_b}$ in the backscattering coefficient is calculated similarly to the error in the scattering coefficient, where the backscattering coefficient $b_b$ is calculated according to Eq. (10), but over the range 90-${150}^{\circ}$.

The relative error of the phase function $p(\theta ) = \beta (\theta )/b$ is calculated as the mean relative error across the range 0.1-${150}^{\circ}$,

The simulated phase functions are measured at log-spaced angles for $\theta _i< {15}^{\circ}$, similar to LISST-VSF measurements. To avoid bias in the forward scattering direction, the data was interpolated to evenly spaced angles with a step size of ${1}^{\circ}$, so that $K=150$.

The errors are analyzed in terms of the optical depth $\tau =cL$ and asymmetry parameter $g$. Here, $L={15}\;\textrm{cm}$ is the length of the cylindrical sample volume. As all samples are simulated with absorption $a={0}\;\textrm{m}^{-1}$, we have $\tau =cL=bL$. The asymmetry parameter $g$ is a measure of how forward peaked the phase function is, and is calculated as the mean cosine of the phase function [29]

The asymmetry parameter ranges from $g=0$ for isotropic scattering and approaches $g=1$ for very forward scattering phase functions.

3. Results and discussion

3.1 Efficiency of VRMs

The efficiency of the VRMs is tested by running the simulation for ${1000}\;\textrm{s}$, both with and without the VRMs. A FF phase function with parameters $n=1.05$ and $m=3.5$ was used for this purpose, simulated with scattering coefficients $b=0.05$, 0.17, 0.55, 1.8, 6, and ${20}\;\textrm{m}^{-1}$. This phase function is used due to its large asymmetry parameter (low scattering in the backward direction), so that the increase in efficiency is expected to be large. In order to evaluate the efficiency of the different VRMs, the VSF obtained with the different VRMs are compared to a reference VSF, which is obtained using no VRMs and a very long simulation time (t>24 hours). The resulting simulated VSFs from the efficiency test are presented in Fig. 3. The parameters used for the simulations with VRMs are $\epsilon =0.2$ and $w_{\textrm {max}}=1$.

 

Fig. 3. Simulated VSF for a total simulation time of 1000s (black) compared to a reference VSF which is simulated without any VRMs for a very long time (t>24 hours). The VSF is plotted for six scattering coefficients $b$, increasing from 0.05 (bottom) to 20 (top). The insets show the volume scattering function plotted from $0 - {20}^{\circ}$ where both the x- and y-axis are in log scale. (a) Simulation without VRMs. (b) simulation with VRMs.

Download Full Size | PDF

It becomes immediately clear from looking at the plots that the applied VRMs result in major noise reduction. This is especially evident in the backward direction and for low scattering coefficients. Also, the new algorithm (with VRMs) reproduces the results obtained with the old algorithm (without VRMs). As the old algorithm has been validated in a previous study (see [14]), the new algorithm is indirectly validated. The discontinuity in the simulated VSFs seen at $\theta = {15}^{\circ}$ and large scattering coefficients is due to the abrupt decrease in single-scattering path length going from the ring detector ($l={15}\;\textrm{cm}$) to eyeball detector ($l={10.3}\;\textrm{cm}$), where a longer path length results in more multiple scattered photons and larger errors [14].

In order to quantify the improvement, we calculate the average relative error in the VSFs over the entire measurement range as

 

Fig. 4. Average relative error between simulated VSFs with a running time of 1000s and a low-variance simulation used as reference. (a) Comparison of VSFs simulated with no VRMs (blue), MFPTR (yellow), DDIS (purple), DDIS and splitting (green), DDIS and MFPTR (orange), and with all VRMs (black). (b) Zoomed in on the lowest error plots.

Download Full Size | PDF

From Fig. 4(a), one can see that the DDIS is by far the most important method for reducing the variance in the simulated VSFs. In fact, comparing the simulation with DDIS to the simulation without any VRMs and a computation time of 10000s, we see that applying DDIS results in over 10x increased efficiency for the five smallest optical depths. At the smallest optical depths, applying DDIS results in much more than 10x efficiency, while it is similar to 10x efficiency at an optical depth of 3. When no VRMs are applied, the relative error decreases with increasing optical depth. This is due to the increased number of photons interacting with the sample volume, and the increased number of multiple scattered photons which increases the probability of scattering to larger angles, resulting in better statistics in the backward direction. Applying only MFPTR, the relative error is significantly decreased at small optical depths compared to the VSF with no VRMs. This is due to the large portion of photons that would pass straight through the sample volume without interaction. As the optical depth increases, this VRM becomes less efficient.

A cropped version of Fig. 4(a) is presented in Fig. 4(b), making it easier to visually differentiate the plots where DDIS is applied. Comparing splitting to MFPTR (combined with DDIS), one can see that MFPTR is more efficient at small optical depths, approximately halving the relative error for the smallest optical depth compared to when only DDIS is applied. As the optical depth increases, the splitting method becomes more important for the overall efficiency. This is an expected result, as larger optical depths lead to increasing number of scattering events per photon, hence more weight accumulation. Applying all VRMs results in the overall best results, having the lowest or close to lowest relative error for all optical depths. The same values for the parameters $\epsilon$ and $w_{\textrm {max}}$ are used for all optical depths shown in Fig. 4. Thus, the VRMs are not optimized and the efficiency of the different methods is expected to improve by tailoring these parameters to the different optical depths. The phase function is also expected to influence the efficiency of the different VRMs and optimal values of the VRM parameters. For instance, the efficiency improvement is expected to be smaller for phase functions with more scattering in the backward direction, as the simulation without VRMs would perform better.

3.2 Multiple scattering errors

Two examples of simulated phase functions are presented in Fig. 5, representing opposite ends of the spectrum of $g$ for the simulations performed in this study. In both examples, it can be seen that the deviation between simulated and true VSFs increases with increasing scattering coefficient (from bottom spectrum to top spectrum), where the deviation is solely due to multiple scattering.

 

Fig. 5. VSF plotted for simulation with all VRMs applied (black) and true values (green). The VSF is plotted for six scattering coefficients $b$, increasing from 0.05 (bottom) to 20 (top). (a) FF phase function with parameters $n=1.05$, $m=3.5$ (b) FF phase function with parameters $n=1.30$ and $m=4.5$. The insets show the volume scattering function plotted from $0- {20}^{\circ}$ where both the x- and y-axis are in log scale.

Download Full Size | PDF

3.2.1 Phase function dependency

In order to analyze how the multiple scattering error depends on the optical depth and phase function, the errors $e_b$ and $e_{b_b}$ and $e_{p}$ are plotted as a function of $\tau$ and $g$. A contour plot for $e_b$ is presented in Fig. 6, where all FF phase functions investigated are included. A second x-axis is added to the top of the graph, displaying the x-values in terms of the scattering coefficient $b$. The phase functions simulated have parameters $m$ ranging from 3.5 to 4.5 in steps of 0.5, and parameters $n$ ranging from 1.05 to 1.35 in steps of 0.05, resulting in a total of 21 phase functions. The phase functions are simulated with the same scattering coefficients as for the VRMs test.

 

Fig. 6. Contour plot of the relative error $e_b$ in the scattering coefficient. The error is plotted as a function of optical depth $\tau$ (or scattering coefficient $b$) and asymmetry parameter $g$. The x- and z-axis (color-axis) are plotted in log scale, while the y-axis is linear. The simulation data is marked as white dots. The black lines indicate a 1%, 10% and 100% relative error, respectively, going from left to right.

Download Full Size | PDF

From Fig. 6, one can see how the relative error increases diagonally from small $\tau$ and $g$ to large $\tau$ and $g$. As expected, the multiple scattering error increases with increasing $\tau$, but there is also a considerable increase with $g$ as well. For the largest $g$, corresponding to $n=1.05$ and $m=3.5$, the relative error in $b$ increases from 1.4% at $\tau =0.0075$ ($b=0.05$) to 494% at $\tau =3$ ($b=20$). An error of 10% is observed at approximately $\tau =0.15$ ($b=1$) while the error is above 100% for $\tau >1$ ($b>7$). Looking at the dependency on $g$ we can see that, for $\tau =3$, the relative error increases from 55% at $g=0.34$ to 494% at $g=0.96$. The dependency on $g$ can largely be explained by the LISST-VSF geometry [14]. Because the radius of the cylindrical sample volume is only ${2}\;\textrm{cm}$, photons scattered at large angles may exit the sample volume with relative short path lengths. On the other hand, photons scattered at small angles may travel along the z-axis, which has a length of ${15}\;\textrm{cm}$. Thus, very forward peaked phase functions will have more photons traveling a longer path length within the sample volume, resulting in more multiple scattered photons and larger errors.

Some distortions can be observed in the top-left corner, caused by a sudden increase in the relative error for some of the phase functions. This is especially apparent at $g=0.83$ and $g=0.91$, where sharp horizontal lines can be seen. Here, there are two different phase functions but with a similar $g$, such that they cannot be visually separated in the plot. It is not clear why this abrupt change in relative error only occurs for small optical depths. However, on closer inspection we find that all the phase functions with increased relative error have the smallest parameter $m=3.5$, where a small $m$ is related to a large VSF in the extreme forward direction. This suggests that extremely forward peaked phase functions are more sensitive to multiple scattering errors at small optical depths. Similarly to Fig. 6, contour plots are plotted in Fig. 7 and 8 for $e_{b_b}$ and $e_p$, respectively.

 

Fig. 7. Contour plot of the relative error $e_{b_b}$ in the backscattering coefficient. The error is plotted as a function of optical depth $\tau$ (or scattering coefficient $b$) and asymmetry parameter $g$. The x- and z-axis (color-axis) are plotted in log scale, while the y-axis is linear. The simulation data is marked as white dots. The black lines indicate a 1%, 10% and 100% relative error, respectively, going from left to right.

Download Full Size | PDF

 

Fig. 8. Contour plot of the relative error $e_p$ in the phase function. The error is plotted as a function of optical depth $\tau$ (or scattering coefficient $b$) and asymmetry parameter $g$. The x- and z-axis (color-axis) are plotted in log scale, while the y-axis is linear. The simulation data is marked as white dots. The black lines indicate a 1% and 10% relative error, respectively, going from left to right.

Download Full Size | PDF

The plot of $e_{b_b}$ displays the same trend seen for $e_b$ in Fig. 6, with $e_{b_b}$ increasing diagonally from small $\tau$ and $g$ to large $\tau$ and $g$. For the largest $g$, the relative error increases from 0.6% to 828%, while for the largest optical depth, the error increases from 74% to 828%. Generally, the error $e_{b_b}$ is larger than $e_b$, which can be explained by the large difference in phase function in the forward and backward direction. As more photons are scattered multiple times, the angular distribution becomes more random. Thus, multiple scattered photons are more evenly distributed over the angular spectrum than the single scattered photons. Adding a photon to an angle where the single-scattering signal is low results in a larger relative increase in signal than adding a photon to an angle where the single-scattering signal is high. Thus, the relative increase in detected photons is larger in the backwards direction, resulting in larger backscattering errors. The distortions seen in Fig. 6 are not observed for the backscattering coefficient, indicating that the distortions are indeed related to forward scattering.

The same diagonal trend is observed for the error $e_{p(\theta )}$, as for $e_b$ and $e_{b_b}$. In general, the error $e_{p}$ is much smaller than the errors in $e_b$ and $e_{b_b}$. This is because the phase functions is normalized, such that the error is only related to the shape of the VSF and not the magnitude. For the largest $g$, the relative error increases from 1% to 45%, while for the largest optical depth, the error increases from 9% to 45%. For $\tau <0.1$, the relative error $e_{p}$ does not increase much, which is due to noise being the dominating source of error when the difference between simulated and true VSF becomes small. This is not the case for $e_b$ and $e_{b_b}$, as both $b$ and $b_b$ are calculated as integrals over an angular range. Thus, fluctuations in the signal are evened out. The plot in Fig. 8 also displays the same distortions seen in Fig. 6, most notably at $g=0.91$.

 

Fig. 9. Relative error $e_b$ in the scattering coefficient plotted as a function of asymmetry parameter $g$. The error is plotted for a variety of FF phase functions with parameter $m=3.5$ (blue), $m=4.0$ (red) and $m=4.5$ (yellow), in addition to a set of Henyey-Greenstein functions (purple).

Download Full Size | PDF

In order to better demonstrate the effect of the phase function on the multiple scattering error, the relative error $e_b$ is plotted in Fig. 9 as a function of only $g$ at the largest optical depth $\tau =3$. Here, we have also included a set of Henyey-Greenstein phase functions, which are much less peaked in the extreme forward direction [29].

From the plot, we can see that the different FF phase functions follow slightly different trend lines, depending on the parameter $m$. For the phase functions investigated here, the phase functions with $m=4.0$ have the smallest error for similar asymmetry parameters $g$, while phase functions with $m=3.5$ and $m=4.5$ seem to have similar errors. However, this cannot be concluded due to the small overlap in $g$ for the $m=3.5$ and $m=4.5$ phase functions. Comparing the FF to the Henyey-Greenstein phase functions, we see that the errors are significantly smaller for the latter. This is due to the very forward peaked nature of the FF phase functions, compared to Henyey-Greenstein phase functions. The difference in error between Henyey-Greenstein and FF phase functions is largest for small $g$, and decreases with increasing $g$. The results presented here demonstrates that, while the asymmetry parameter $g$ is a good indicator for the size of the error, details in the phase function may influence the error significantly.

3.2.2 Absorption dependency

The results presented in the previous section were all obtained with an absorption coefficient of $a={0}\;\textrm{m}^{-1}$. In Fig. 10, the error $e_b$ is plotted for a phase function simulated with $a=0c$, $a=0.2c$ and $a=0.4c$, in order to investigate the effect of absorption. The same set of attenuation coefficients, i.e. optical depths, is used for each absorption coefficient. Thus, the only difference between each simulation set is the ratio of scattering to absorption.

 

Fig. 10. The relative error $e_b$ in the scattering coefficient plotted as function of (a) attenuation coefficient $c$, and (b) scattering coefficient $b$. The error is plotted for an absorption of $a=0c$ (blue), $a=0.2c$ (orange) and $a=0.4c$ (yellow).

Download Full Size | PDF

In Fig. 10(a), the error in the scattering coefficient $e_b$ is plotted against the attenuation coefficient $c$. Here, we see that $e_b$ decreases with increasing absorption, which is expected as the scattering decreases, hence the number of multiple scattered photons decreases. The difference in error between the different values of $a$ increase with increasing optical depths. At the smallest optical depth, the errors at $a=0c$ and $a=0.4c$ are $1.43\%$ and $1.28\%$, respectively, while at the largest optical depth, the errors are $494\%$ and $180\%$. In Fig. 10, $e_b$ is plotted against the scattering coefficient $b$. Here, it is clearly demonstrated that the multiple scattering error is only dependent on $b$. Thus, when evaluating the multiple scattering error in a measurement, it is better to use the scattering coefficient, or the size $bL$, rather than the attenuation coefficient or the optical depth $\tau =cL$.

4. Summary and conclusion

We have performed an extensive investigation into the errors originating from multiple scattering for the LISST-VSF instrument. Several different Henyey-Greenstein and FF phase functions have been simulated, with scattering coefficients $b$ ranging from 0.05 to ${20}\;\textrm{m}^{-1}$, and with asymmetry parameters $g$ ranging from 0.34 to 0.96. The results show that errors in the scattering coefficient $b$, backscattering coefficient $b_b$, and phase function $p(\theta )$, depend heavily on both the scattering coefficient and asymmetry parameter. We find that multiple scattering errors may only be considered negligible (<$10\%$) when the scattering coefficient is smaller than 1 m$^{-1}$. The largest error is observed for the backscattering coefficient, where the error ranges from 74% to 828% for the different phase functions at an optical depth of $\tau =3$. For the scattering coefficient, the error ranges from 55% to 494%, and the phase function from 9% to 45%. The results also show significant differences between Henyey-Greenstein and FF phase functions, where the error is considerable larger for the FF phase functions for similar asymmetry parameters. We conclude that this is due to the extremely forward peaked shape of the FF functions compared to the Henyey-Greenstein functions, resulting in a longer average distance traveled within the sample volume and more multiple scattered photons. This observation demonstrates the importance of the shape of the phase function when measuring the VSF with the LISST-VSF instrument. The effect of absorption was also investigated. We found, when the error is plotted against optical depth, the error decreases with increasing absorption. However, when plotted against the scattering coefficient, the error was independent of absorption, i.e. the error is only dependent on the scattering coefficient for a given phase function.

The LISST-VSF has an angle-dependent path length varying between 10.3 and 17.5 cm, which accentuates multiple scattering effects on the measured VSF in turbid waters, in particular after applying the absolute calibration [13]. This does not mean other VSF meters are not affected by multiple scattering. The MASCOT and MVSM instrument, two VSF meters that have been widely used for modern era VSF measurements, have reported constant path lengths of 20.0 cm [30,31]. The longer path lengths imply that the multiple scattering errors can be significant for $b>1$ m$^{-1}$. This can have implications for historical data sets in waters with moderate-to-high turbidity. However, we should note that the results in this study are instrument-specific to the LISST-VSF.

A set of variance reduction methods has been applied to the previously developed Monte Carlo simulation, resulting in drastically reduced computation time. This enables a fast method for analysing errors in the measured VSFs. The most important of these methods is detector directional importance sampling (DDIS), where the photons have an increased probability of being scattered towards the eyeball detector. The applied methods result in over 10 times increased efficiency for optical depths below 0.9, and similar to 10 times increased efficiency up to an optical depth of $\tau =3$. The efficiency is expected to be further enhanced by optimizing the parameters used for the variance reduction methods. The efficiency is also expected to depend on the phase function, where the improvement in efficiency is expected to increase with increasing asymmetry parameter. This is because the simulation without variance reduction methods achieves a very low signal in the backward direction when the asymmetry parameter is large. For small asymmetry parameters, many photons are naturally scattered in the backward direction, reducing the need for DDIS. The implementation of the detector directional importance sampling in this simulation is geometry specific, however, the method of implementation might be relevant for other cases.

For our future work we aim to correct the errors arising from multiple scattering in the LISST-VSF instrument. The main goal is to correct the VSF over the entire LISST-VSF measuring range, thus correcting both the scattering coefficient, backscattering coefficient, and phase function at the same time.