2.1 Vertical profiles

The molecular and aerosol profiles are described below. It should be noted that we assume the atmosphere to be plane-parallel; i.e. there is no correction applied for sphericity of the Earth's atmosphere.

2.2 Aerosol types

Four aerosol types have been defined. Aerosol size distributions are assumed to be log-normal with the parameters shown in Table 1. Refractive indices are taken from the literature (see Table 1) and the aerosols are assumed to be internally mixed, with a volume-weighting rule. The four aerosol types defined in this benchmark have been chosen to represent a variety of size and optical properties observed in the atmosphere. Two fine-mode aerosols are defined: industrial scattering (aer_is) is sulfate aerosol with a relative humidity of 50 %, while industrial absorbing (aer_ia) is a mixture of the Optical Properties of Aerosols and Clouds (OPAC; Hess et al., 1998) aerosol type INSO (insoluble), SOOT (soot), and WASO (water soluble, at 50 % of relative humidity). The industrial absorbing aerosol is similar to an organic-matter aerosol type. We also consider a monomodal dust aerosol and a bimodal sea salt aerosol (aer_ss).

Table 1Aerosol physical and optical properties for the different aerosol types considered in the benchmark. σ_g is the geometric standard deviation. The refractive indices are linearly interpolated from the OPAC database (Hess et al., 1998). Dust refractive indices are linearly interpolated from Woodward (2001).

Download Print Version | Download XLSX

2.3 Aerosol optical properties

Aerosol optical properties were computed by the authors using their own Mie routines, and results were cross-checked for consistency. The Mie routine used to generate the input data was rewritten in Fortran based on the implementation of Toon and Ackerman (1981) and was compiled with the Intel Fortran compiler in quadruple precision. The calculations are monochromatic (i.e. no integration was performed on the wavelength). The aerosol size distribution is integrated over the size range reported for each type in Table 1, with a Gaussian quadrature of 10 000 points in logarithm of the radius. The number of terms used to compute the sums of the a_n and b_n coefficients is approximated as a function of the Mie parameter ($x=\mathrm{2}\mathit{\pi }r/\mathit{\lambda }$) following Wiscombe (1980).

For each aerosol type, we provide the aerosol optical depths at the 440, 670, 865, and 1024 nm wavelengths corresponding to the selected AOD at 550 nm. The aerosol mass extinction coefficient (although not necessary to conduct the RTM calculations), single-scattering albedo, and asymmetry parameter are also given.

The phase matrix is computed and provided at 50 000 points evenly spaced in scattering angle from 0 to 180^∘. As the particles are assumed to be spherical, only the S₁₁, S₁₂, S₃₃, and S₃₄ elements of the phase matrix are non-zero and provided. We also decompose the phase matrix in generalized spherical functions (de Rooij and van der Stap, 1984; de Haan et al., 1987) as needed for the vector version of the discrete ordinate method. For this, the phase matrix is integrated with an 8000-point Gaussian quadrature over the cosine of the scattering angle, and we provide the first 750 moments of that decomposition for the six commonly used elements (α₁, α₂, α₃, α₄, β₁, β₂). For completeness, we also provide the corresponding Legendre coefficients (by definition, the lth Legendre coefficient equals the lth Legendre moment multiplied by 2l+1). Thus, we provide all necessary information on aerosol optical properties for the usual RTMs so that no further processing of the aerosol optical properties should be required. Of course, the RTMs that neglect polarization should only consider the S₁₁ term of the phase matrix or the α₁ Legendre moments.

The accurate reconstruction of the original phase function from the Legendre decomposition (S₁₁) cannot be ensured when a low number of Legendre moments is used to approximate it. As for the vertical resolution, there is no obligation for this benchmarking exercise to represent the exact phase function, and it may be interesting to approximate it to strike a balance between accuracy and computational cost (see, for example, Ding et al., 2009). The minimum number of Legendre moments for which the approximated phase function is positive depends on the aerosol type and the wavelength (basically, on how large the forward peak is). For the industrial scattering and industrial absorbing aerosols, a minimum of four to eight moments is required (depending on the wavelength). For dust aerosol, the minimum is between 14 moments (at 1024 nm) and 40 moments (at 470 nm), while for sea salt, the minimum is between 144 moments (at 1024 nm) and 348 moments (at 470 nm).

Figure 1 shows the phase function and the linear depolarization ratio for the four aerosol types considered in this benchmark.

Figure 1First (i.e. S₁₁, left panels) and second (i.e. S₁₂, middle panels) element of the phase matrix and linear depolarization ratio (i.e. $-S_{\mathrm{12}}/S_{\mathrm{11}}$) (right panels) as a function of the scattering angle (^∘) for the four aerosol types (rows), where aer_du stands for dust, aer_ia for industrial absorbing, aer_ss for sea salt, and aer_is for industrial scattering. The colour code corresponds to the five wavelengths (470, 550, 660, 865, and 1024 nm). The mass extinction coefficient (ext, in m² g⁻¹), single-scattering albedo (ssa), and asymmetry parameter (g) are also indicated in each panel. The phase function is normalized to 2 $\left(\text{i.e.}{\int }_{\mathrm{0}}^{\mathit{\pi }}{S}_{\mathrm{11}}\left(\mathrm{\Theta }\right)\sin \mathrm{\Theta }\mathrm{d}\mathrm{\Theta }=\mathrm{2}\right)$.

Download

2.4 Surface reflectance

Two surface reflectance models are selected: a Lambertian and an oceanic bidirectional reflectance distribution function (BRDF) model. For the Lambertian model, three surface albedos were chosen, namely 0, 0.05, and 0.1. The spectral dependence of the surface reflectance is ignored; i.e. the albedo is assumed to be the same for the five wavelengths considered. The second surface reflectance model is the oceanic glitter model from Mishchenko and Travis (1997). This model is suited for both full vector radiative transfer models (i.e. including polarization) and for the scalar approximations.

Wind speed is set to 10 m s⁻¹, but the Mishchenko and Travis (1997) code assumes an isotropic Gaussian distribution of (oceanic) surface slopes, so the wind orientation is not taken into account. The salinity is set to 34.3 ‰ and water refractive indices are interpolated, for each wavelength, from the data provided in Hale and Querry (1973) as inputs to the Mishchenko and Travis (1997) routine (available at https://www.giss.nasa.gov/staff/mmishchenko/brf/, last access: 19 February 2019). The effect of shadowing is also taken into account.

The convention for the azimuthal angle used in the Mishchenko and Travis (1997) routine is the opposite of the convention that we use below to define geometries, so the coupling between the atmospheric RTM and the oceanic code has to be done carefully.

2.5 Geometries and configurations

We have selected eight cases for aerosol loadings in the atmosphere, with AOD at 550 nm ranging between 0 (no aerosol) and 2 (high aerosol loading). These values refer to the 550 nm wavelength, and the corresponding AODs can be estimated for the other wavelengths depending on the aerosol types. A table with the equivalences in AOD for each aerosol type and wavelength is also provided in the benchmark protocol.

We have further selected four solar and four viewing zenith angles – θ_s and θ_v – from 0^∘ (zenith) up to 60^∘. Given the provided surface reflectance models and the assumption of spherical (hence randomly oriented) aerosols in the atmosphere, the radiance at the top of the atmosphere is only dependent on the difference between the incident and reflected azimuthal angle and not on each value taken independently. We have chosen five azimuthal angles between ϕ_r=0^∘ (backscatter when solar and viewing zenith angles are equal) and ϕ_r=180^∘ (forward scattering when solar and viewing zenith angles are equal). The different cases are summarized in Table 2. Removing the redundancies in aerosol types for AOD =0 and azimuthal angle for ${\mathit{\theta }}_{\mathrm{s}}={\mathit{\theta }}_{\mathrm{v}}={\mathrm{0}}^{\circ }$, this results in 11 020 cases for the BRDF surface and 33 060 cases for the Lambertian surfaces (see Table 3), which gives a grand total of 44 080 cases.

Table 2Summary of the proposed geometrical and environmental conditions used in the simulations.

Download Print Version | Download XLSX

Table 3Statistical measures of DISORT and scalar VLIDORT reflectances against vector VLIDORT (first row). Errors of the last rows are defined as the absolute value of the relative error, using the VLIDORT (vector) reflectance as the reference. VLIDORT is the configuration with 32 streams and the TMS correction activated. N refers to the number of cases.

Download Print Version | Download XLSX

2.6 Measured variables

We aim to provide a standardized dataset to evaluate the possibility of assimilating satellite radiances in operational forecasts of atmospheric composition. Thus, two variables are important to estimate for each case: the simulated reflectance (or radiance) that would be observed by the satellite and a “sound” estimate of the radiative transfer code runtime.

We will describe the accuracy of the models in terms of the reflectance, ρ_I, of the intensity (i.e. the first Stoke component), which is related to the radiance, I, through the following relationship:

where I_s is the incident solar flux at the top of the atmosphere (conveniently taken to equal 1 for reflectance calculations) and θ_s is the solar zenith angle.

The runtime may be difficult to estimate and compare between different models. We aim to measure the runtime of the computations excluding input/output operations, as typical applications would embed the RTM code in a larger operational system. The runtime of the RTM should be measured in a single-thread configuration. Ideally, the simulations should be performed on the same core of the computer and with the same system load conditions (except for the RTM) as other threads on the operating system (in particular with a high input/output load) could hamper the correct timing of the model. Multiple repetitions of the computations can be performed to get better estimates.

Given a set of atmospheric and aerosol conditions, some models can compute, almost without extra cost, the radiances for several viewing geometries, while other models require multiple simulations to achieve the same. For this protocol, this would introduce a difference of a factor of circa 80 (the number of viewing geometries) in the reported runtimes. Depending on the application, this feature of some models may or may not be relevant. For this reason, we present here the runtimes of the tested models in two configurations: one taking advantage of the simultaneous output for multiple viewing geometries and another one without this feature. In the second configuration (called “ind” in Sect. 3), each simulation outputs only one geometry while in the first configuration (called “mult” in Sect. 3), it outputs as many values as possible within the specifications of the protocol. The number of multiple outputs per simulation will be detailed in the next section.

3.1 Reference model: VLIDORT

The VLIDORT model is an independent implementation of the discrete ordinate method based in part on earlier work by Siewert (2000a, b) and other sources. In this work we used the VLIDORT version 2.7, which includes a BRDF supplement that can be called before the main program. A detailed user's guide is available (Spurr, 2014). One interesting feature of this model is that it also provides Jacobians with respect to the model inputs (not used in this study). In this model, it is possible to choose how many elements are computed in the Stokes vector. Although VLIDORT is able to solve the radiative transfer equation for the four-element Stokes vector $(I,Q,U,V)$ (Spurr, 2008), we will refer to the vector runs when the 3×3 approximation of the RT problem is used (i.e. neglecting circular polarization in the atmosphere by solving only for $(I,Q,U)$) and we will refer to the scalar runs when only the first element of the Stokes vector is used in the computations (i.e. when polarization of electromagnetic radiation is not accounted for). VLIDORT is written in double precision; it has the capability to output, in a single call, radiances for any number of solar and viewing angles. The protocol vertical resolution with 49 layers is used. It should be noted that the discrete ordinate method shows a numerical instability for zenith viewing angle θ_v of 0^∘ or a relative azimuthal angle ϕ_r of 90^∘. These angles are therefore prescribed at θ_v=0.0001^∘ and ϕ_r=89.98^∘ instead to circumvent this problem. We use the model with 32 streams, and the TMS correction (Nakajima and Tanaka, 1988) is activated. The TMS correction removes the single-scattering feature computed with a few streams and replaces it by the exact single-scattering contribution. This correction helps to remove non-physical features due to the use of a low number of Legendre coefficients in the computation. The VLIDORT vector is our reference model.

3.2 DISORT

The DISORT model (Stamnes et al., 1988) is a flexible and widely used RTM. It solves the scalar approximation of the radiative transfer equation in a plane-parallel atmosphere by using the discrete ordinate method. In this study, we use the DISORT code, version 4, available from http://lllab.phy.stevens.edu/disort/ (last access: 19 February 2019), and we coupled the oceanic BRDF following the code developer instructions. We use the model with 32 streams, and the TMS correction (Nakajima and Tanaka, 1988) is activated. It should be noted that most of this code is written in single precision, but important routines (such as the computation of the Gaussian nodes) are written in double precision. We use 32 streams in our main configuration, but sensitivity results for 4, 8, 16, 32, 64, 96, and 128 streams are also presented below. The protocol vertical resolution with 49 layers is used. Similar to VLIDORT, the DISORT code has numerical instabilities which manifest themselves when the solar and viewing zenith angles (θ_s and θ_v) are equal. In such cases, we perturb the cosine of the viewing zenith angle by $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{5}}$.

Given a solar zenith angle, the DISORT model can provide outputs for multiple viewing geometries. In this study, the number of geometries for the mult timing configuration is the number of ϕ_r cases multiplied by the number of θ_v cases (i.e. $\mathrm{4}\times \mathrm{5}=\mathrm{20}$).

Although we know the DISORT model to be too slow for online data assimilation of aerosol radiances, we find it useful to document here the trade-off between accuracy and computing time as a function of the number of streams used in the calculations. Figure 2 shows both the mean fractional error and the average computing time for 4, 8, 16, 32, 64, 96, and 128 streams. The mean fractional error is shown relative to our reference VLIDORT (vector) and is of the order of 1 or 1.5 % for the Lambertian and oceanic BRDF surface cases, respectively, as soon as the number of streams exceeds eight. Most of that discrepancy is because DISORT neglects the polarization. When DISORT is compared to VLIDORT (scalar), the mean fractional error is significantly smaller: <0.1 % for the Lambertian surface case and <1 % for the surface BRDF case.

There is a big gain in accuracy when increasing the number of streams from 4 to 8 and little further gain in accuracy beyond 16 streams. A deterioration in accuracy is observed beyond 96 streams in the case of the BRDF surface, which may be due to the DISORT model being coded in single-precision float. While the dependence of model accuracy on the number of streams is well known (Stamnes and Swanson, 1981), we find it useful to confront the loss of accuracy with the gain in computational cost. Computing time increases rapidly with the number of streams (faster than the cube of the number of streams beyond 16 streams). Using 16 streams might, therefore, be an interesting trade-off between accuracy and computational cost.

Figure 2Accuracy (in black) and computing times (in red) for the DISORT model as a function of the number of streams used. The Lambertian and oceanic BRDF surface cases are shown with solid and dashed lines, respectively. The accuracy is shown in terms of mean fractional error against VLIDORT (vector). The computing times are an average for 20 geometries and were estimated on a processor AMD Opteron 6378, 2.4 GHz. Please note the logarithmic scales for the number of streams used in the DISORT calculations (x axis) and the timings (right vertical axis).

Download

3.3 6SV2

The 6SV2 model is a radiative transfer model based on the successive order of scattering (SOS) method (Vermote et al., 1997). It solves the radiative transfer equation by using the three-element approximation (3×3) of the Stokes vector. The accuracy of the code against previous benchmark and Monte Carlo codes has been reported to be better than 1 % (Kotchenova et al., 2006; Kotchenova and Vermote, 2007). We have used version 2.1 of the model. Two configurations are available: a standard accuracy configuration that consists of 83 Gaussian points for the phase function, 30 vertical levels, 25 zenith angles, and 181 azimuth angles and a high-accuracy configuration that consists of 121 Gaussian points for the phase function, 50 vertical levels, 75 zenith angles, and 181 azimuth angles. We will present results only from the standard accuracy, which is 2 to 70 faster than the high-accuracy configuration, but will comment briefly on the high-accuracy configuration in Sect. 4.2.

Molecular scattering and absorption vertical profiles are embedded in the code. Both are internally computed according to the atmospheric profile as part of “core” SOS routines. It is therefore not appropriate to prescribe our own molecular scattering and absorption profiles, so we use the 6SV2 defaults instead. The difference in the total molecular scattering optical depth is less than 0.26 % with respect to the files provided in this protocol. We have coupled the BRDF calculations to the 6SV2 model, but it should be noted that the atmosphere–surface coupling of the oceanic BRDF model in 6SV2 does not use the full surface reflection matrix but only the first column of it (i.e. the column that relates reflected I with incoming I, Q, and U).

The 6SV2 model can only compute one viewing geometry in each model call (i.e. the ind configuration is exactly the same as the mult configuration).

3.4 FLOTSAM

We also test in this study a preliminary version of the Forward-Lobe Two-Stream Radiance Model (FLOTSAM), which is being developed at the European Centre for Medium-Range Weather Forecasts (ECMWF). Unlike a lookup-table-based solar RTM, FLOTSAM permits an arbitrary layering of different types of scatterers, including clouds, aerosols, and molecules. The model has been designed to be fast enough to be used in iterative assimilation and retrieval schemes. It exploits the fact that particle phase functions typically contain a “narrow forward lobe” associated with scattering angles of less than a few degrees and a “wide forward lobe” associated with scattering angles of the order of 15^∘. This enables coupled ordinary differential equations to be written down for radiation in the quasi-direct solar beam (including the narrow forward lobe) and a wider beam of radiation propagating away from the Sun due to scattering by the wide forward lobe. Additionally, the upwelling and downwelling diffuse fluxes are computed using the two-stream approximation. FLOTSAM may, therefore, be thought of as the solar analogue of the infrared two-stream source function technique of Toon et al. (1989). To facilitate the use of FLOTSAM in variational assimilation and retrieval schemes, it has been coded in C++ using version 2.0 of the Adept library (Hogan, 2014, 2017), which enables the Jacobian to be easily computed via Adept's fast automatic differentiation capability. FLOTSAM has been incorporated into the Cloud, Aerosol and Precipitation from Multiple Instruments using a Variational Technique (CAPTIVATE) retrieval scheme, which will provide operational products from the forthcoming EarthCARE satellite (Illingworth et al., 2015).

FLOTSAM is designed to be used in iterative schemes in which multiple calculations are performed with different profiles of optical properties but for the same viewing geometry. This is implemented by allowing the user to perform the set-up calculations once for a particular solar zenith angle, viewing zenith angle, and viewing azimuth angle and then to reuse the information in subsequent radiance calculations for different optical profiles. We assess the computational cost for both the ind and the mult configurations, with the latter consisting of $\mathrm{4}\times \mathrm{4}\times \mathrm{5}\times \mathrm{8}=\mathrm{640}$ cases involving 80 geometry set-up calculations and each separate geometry being used with eight different optical depths.

4.1 Methods and statistics

The protocol has been tested for the five model versions described above: VLIDORT (vector), VLIDORT (scalar), DISORT, 6SV2, and FLOTSAM. As explained in Sect. 2.6, the accuracy is quantified in terms of monochromatic reflectances at the top of the atmosphere and the accuracy is computed against the vector version of the VLIDORT model unless stated otherwise.

With as many as 44 080 cases considered in this study, we have to find original ways to visualize the results as explained below. Along with the average reflectance computed for each model, the following statistical measures have also been computed:

where N is the number of experiments, m_i is the reflectance of the tested model, and r_i is the reflectance of the reference model for case i. By definition, the mean fractional bias can range between −2 and 2, with a best score of zero, while the mean fractional errors can range between 0 and 2, with a best score of zero. Both values are multiplied by 100 with respect to the definitions above in the tables presented below. We also state the percentages of cases where the relative error of the reflectances is above defined relative error thresholds. These diagnostics, which could be detailed for particular subsets of our 44 080 cases, can help to decide if and under which conditions a particular RTM could be useful, according to the required accuracy of the reflectance calculations.

Runtimes of the models are presented in an absolute measure (seconds), and no reference is needed. Finally, we will comment on possible trade-offs between accuracy and computing time.

4.2 Accuracy

Figure 3 shows histograms of the simulated reflectances of the VLIDORT (vector) model at the top of the atmosphere for all or various subsets of the 44 080 cases. These reflectances form the reference simulations for this study. The width of the histograms (or violins) in Fig. 3 is proportional to the number of cases in each bin of the reflectance (y axis) discretization. The multiple panels (and multiple violins) are computed with all the simulations but filtered by the value indicated on the x axis. Therefore, they can be interpreted as analogous to marginal density probability functions of the set of simulations. Figure 3 shows a wide range of reflectances (from near 0 to almost 0.5) and known dependencies of the reflectance with respect to some of the input variables such as the surface albedo (Lambertian case), the AOD, or the wavelength in this range of the spectrum.

Figure 3Histograms of the computed reflectances for the VLIDORT model. Panel (b) is a histogram of all considered cases, while the other panels classify the cases according to the value of θ_s (a), θ_v (c), ${\mathit{\varphi }}_{\mathrm{r}}={\mathit{\varphi }}_{\mathrm{v}}-{\mathit{\varphi }}_{\mathrm{s}}$ (d, “no ϕ_r” corresponds to the cases where ${\mathit{\theta }}_{\mathrm{s}}={\mathit{\theta }}_{\mathrm{v}}=\mathrm{0.0}$), and surface type – BRDF or Lambertian (e), surface albedo value (f), 550 nm aerosol optical depth (g), wavelength (in µm, h), and aerosol type (“no aer” corresponds to AOD =0, i). The width of the classes is indicated at the bottom right ($\mathrm{\Delta }y=\mathrm{2.49}\times {\mathrm{10}}^{-\mathrm{3}}$), and the maximum number of cases in a class is provided for each histogram. The histograms are presented as “violin plots” with a rotated density plot on each side. Mean (triangles) and median (circles) values are indicated. Note that the scale used on the x axis is different for each violin plot.

Download

Figure 4Same as Fig. 3 but for the relative error of the DISORT and VLIDORT (scalar) models against the vector version of VLIDORT. Please note the ±5 % error range on the y axis. The number of cases outside the range is indicated in Table 3. The width of the classes is $\mathrm{\Delta }y=\mathrm{4.98}\times {\mathrm{10}}^{-\mathrm{4}}$. Note that the scale used on the x axis is different for each violin plot but is the same for the two models within a violin plot.

Download

A first comparison of the relative error is shown in Fig. 4 for two scalar models (with otherwise the same degree of sophistication as the reference model). Figure 4 shows the histograms of the relative errors, $(m-r)/r$, with m being the tested model and r the reference model (the vector version of VLIDORT). It can be seen that for both models, DISORT and the scalar version of VLIDORT, the errors are mostly below 5 %. They are centred on zero for the Lambertian surface, but a negative bias is observed for the BRDF case. Comparison of scalar VLIDORT against vector VLIDORT provides a direct measure of the relative errors due to neglecting polarization in the model. As expected and because of a larger molecular scattering contribution, we observe larger relative errors for the shorter wavelengths (MFE =1.8 % at 470 nm; MFE =0.5 % at 1024 nm). We also observe larger errors for the accumulation-mode aerosols (industrial scattering and industrial absorption have MFE ≈1.2 %) in contrast to coarse-mode aerosols (dust and sea salt have MFE ≈0.6 %). These errors are associated with the larger polarized phase function (see Fig. 1) of industrial scattering and industrial absorption aerosols for large AOD (see Figs. A2 and A4). Although no clear relation is shown for relative errors as a function of the geometrical parameters in Fig. 4, Fig. A2 shows a clear link between the relative errors and the scattering angle for these two accumulation-mode aerosols at large AOD. The oceanic BRDF function used in this study also enhances the errors of the scalar versions as these neglect the polarizing effect of the surface because the glint can increase the polarized feature of the light for some angles. Larger relative errors are shown for the dark surface case (Lambertian surface albedo equal to zero has MFE ≈1.2 %) because the reference reflectances are smaller than for their 0.05 and 0.1 Lambertian counterparts (with MFE ≈0.7 % and MFE ≈0.5 %, respectively). The error statistics of the VLIDORT (scalar) and DISORT models are very similar, with noticeable differences only for the BRDF case (Fig. 4 and Table 3). The scalar approximation may thus be appropriate for some combinations of wavelengths, aerosols types, and surface types but not for others. We further note, that for both types of surfaces, the scalar VLIDORT simulations outperform DISORT simulations when they are compared against the vector version of VLIDORT (which is not surprising), but in both cases the overall mean fractional error is larger than 0.7 % and the mean fractional bias is negative.

Figure 5Same as Fig. 3 but for the relative error of the FLOTSAM and 6SV2 models against the vector version of VLIDORT. Please note the ±10 % relative error range of the y axis. The width of the classes is $\mathrm{\Delta }y=\mathrm{9.95}\times {\mathrm{10}}^{-\mathrm{4}}$. The number of cases outside the range is indicated in Table 4. Note that the scale used on the x axis is different for each violin plot but is the same for the two models within a violin plot.

Download

Figure 6Same as Fig. 3 but for the runtime of the 6SV2 model. Please note the logarithmic scale on the y axis.

Download

Figure 5 and Table 4 show the relative errors of the other two models used in this study, 6SV2 and FLOTSAM, against VLIDORT (vector). As the FLOTSAM model does not include polarization, it is not expected to outperform the scalar approximations exemplified in the DISORT and VLIDORT (scalar) comparisons. However, this is the case for the 470 nm cases, as can be seen from the corresponding violin plot of Fig. 5. For 6SV2, the relative error shows a mode around zero, while for FLOTSAM, the shapes of the violin plots are more similar to those of the scalar approximations shown in Fig. 4. The reason for the unexpected dependence of the 6SV2 model errors on the wavelength is not completely clear, but it could be related to the technical difficulties in prescribing user-specific molecular scattering and absorption profiles in the 6SV2 code (see Sects. 2.1.1 and 3.3). In contrast to Fig. 4, Fig. 5 shows larger relative errors for sea salt aerosols for both models. While larger errors in Fig. 4 were related to the polarization/depolarization effects from aerosol particles (as shown in the right column of Fig. 1), errors of 6SV2 and FLOTSAM could be related to the more pronounced forward-scattering peak of the sea salt aerosols and their large single-scattering albedo. In fact, errors for large AOD (Fig. A3) are not strongly linked to the amount of polarized radiance (Fig. A4) for sea salt aerosols. In the case of 6SV2, the coarse discretization of the aerosol phase function used in the model configuration (83 Gaussian points, Sect. 3.3) hampers the accurate simulation of the sea salt forward-scattering peak. There is no clear dependence of the relative errors on the viewing and solar geometries but, as for Fig. 4, larger errors are observed for larger AOD and for the surface BRDF case.

Table 4Same as Table 3 but for the 6SV2 and FLOTSAM models.

Download Print Version | Download XLSX

Synthetic error statistics are shown in Table 4. Both the 6SV2 and FLOTSAM models have comparable mean fractional errors, but these are larger than those of DISORT and VLIDORT (scalar). Some difference concerning the type of surface can be appreciated, with 6SV2 performing better in the case of Lambertian surfaces. It should be noted that using the high-accuracy configuration of 6SV2 instead of the standard accuracy configuration reduces the mean fractional error from 1.5 % to 1.1 % for the Lambertian surface case but does not improve this score for the BRDF surface case. The mean fractional bias is negative, and it is of similar magnitude for both models. However, 6SV2 presents larger biases over the two types of surfaces (and the sign of the bias differs for the BRDF and Lambertian surfaces). In contrast, FLOTSAM shows virtually no bias for the Lambertian surface cases and a smaller bias than 6SV2 for the BRDF cases. An important difference with the scalar models is the presence of a heavy tail in the relative error distribution, which can be inferred from the last rows of Table 4.

4.3 Timing

In this section, we present statistics of computing times for three of the models used in this study. All the computing times presented below were estimated on an Intel i5-4690, 3.5 GHz workstation running Linux (kernel 2.6.32), GNU C++ and Fortran compiler version 4.4.7, on a single processor and with exclusive use of the workstation. We have compiled the models using their default compilation flags, i.e. “-O” for 6SV2, “-O3” and “-march=native” for FLOTSAM, and “-O3” for DISORT.

We start this section with the 6SV2 model which can only be used in ind configuration. Figure 6 shows the violin plots of the computing time of the 6SV2 model that spans as much as 2 orders of magnitude. Clearly, the BRDF cases are more expensive than the Lambertian surface cases, with most computing times ranging between ca. 0.7 and 4 s for the former. This appears to be due to multiple calls to the BRDF routine. It may be possible to optimize this routine or recode it in a more efficient way, but we have not attempted to do so. We also note an increase in computing time with increasing solar zenith angle (θ_s) and AOD and with decreasing wavelength. This is consistent with a larger number of scattering orders being necessary in such cases. Cases for coarse mode aerosols (especially sea salt aerosols) are more expensive than cases for accumulation-mode aerosols, with some calculations reaching 5 s for sea salt aerosols. We do not have an explanation as to why the cases for industrial absorbing aerosols are more expensive than those for industrial scattering aerosols.

Table 5Average runtimes in seconds. The standard deviation is shown in parentheses. The mult runtimes account for 640 profiles for FLOTSAM, 20 profiles for DISORT, and 1 profile per model call for 6SV2. The one-off set-up costs of FLOTSAM are invoked only once per instance of the executable, no matter how many radiances are subsequently calculated, so are the same for the ind and mult experiments.

Download Print Version | Download XLSX

As discussed above, the DISORT and FLOTSAM models offer efficiency savings when multiple profiles are computed together. In the case of DISORT, multiple viewing geometries can be computed for one solar zenith angle. In the case of FLOTSAM, a smaller saving is that the geometry set-up calculations may be computed once for a given solar/viewing geometry and then reused for multiple atmospheric profiles with different optical properties. The runtime of these models is, in general, not sensitive to the viewing geometry, wavelength, aerosol type, or AOD. However, there are some differences between the two types of surface reflectance. Table 5 shows a summary of the measured runtimes for the three tested models. Computing times are indicated for the two measurement procedures (ind and mult) and for the two types of surface reflectance (BRDF and Lambertian). In the case of FLOTSAM, we also report the “one-off set-up” cost, which in the case of the BRDF involves creating a four-dimensional lookup table. We stress that this is invoked only once per call to the executable, so in an operational data assimilation system performing many millions of radiance calculations it would increase the executable runtime by only a very small fraction. The one-off set-up costs are not shown for the 6SV2 and DISORT models as they are not isolated from the rest of the computations.

For the 6SV2 model, values of Table 5 are in concordance with Fig. 6 with runtimes of the order of seconds. The DISORT model average runtime is smaller than that of 6SV2. Similarly to 6SV2, there is an important dependence of the DISORT runtime on the type of surface considered. The extra cost of computing 20 profiles, instead of 1 profile, can be estimated as the difference between the mult and ind cases. For DISORT, this additional cost is about 0.07 s for the Lambertian case and about 0.22 s for the BRDF case, which correspond to 25 % and 18 % of the total mult runtime. This would benefit an application that requires simultaneous calculation of many viewing geometries. The spread of DISORT runtimes indicated in Table 5 arises because DISORT efficiently decreases the number of internal computations for symmetric geometries (solar zenith angle equal to zero, 21 % of the cases), and it is faster for AOD =0 (3.4 % of the cases). For the remaining cases of Table 5, there is almost no spread in the DISORT runtime.

Table 5 shows the runtimes of the FLOTSAM model. Firstly, it is noted that the standard deviation of the reported values is small (around 1 order of magnitude less than the values). This is because FLOTSAM performs the same number of operations for every case; thus, there is no added value of showing histograms of runtimes. For this model, it can be seen that the one-off set-up cost of the BRDF is a little under 3 s, but in the ind case, subsequent radiance calculations are 3–4 orders of magnitude faster than the other two models tested, and the BRDF calculation is only slightly slower. If many profiles have to be computed with the same observation geometry, then the geometry set-up cost is invoked only once and the average cost per profile is reduced. This is illustrated in the mult case in which 640 profiles are computed, which involves 80 different geometries. Excluding the model set-up time, the computation of a single profile takes around 0.3 to 0.4 ms in the ind case and around $\left(\mathrm{99}\mathrm{ms}\right)/\mathrm{640}=\mathrm{0.16}$ ms in the mult case.