A Functionalized Monte Carlo 3D Radiative Transfer Model: Radiative Effects of Clouds Over Reflecting Surfaces

In the Earth Sciences, the 3D radiative transfer equation is often solved for by Monte Carlo (MC) methods. They can, however, be computationally taxing, and that can narrow their range of application and limit their use in explorations of model parameter spaces. A novel family of MC algorithms is investigated here in which single simulations provide estimates of both radiative quantities A for a set of parameters x^ $\hat{x}$ , as usual, as well as the overarching functional A $\mathcal{A}$ (x) that can be evaluated, extremely efficiently, at any x. One such algorithm is developed and demonstrated for horizontally averaged broadband solar radiative fluxes as functions of surface albedo for uniform Lambertian surfaces beneath inhomogeneous cloudy atmospheres. Simulations for a high‐resolution synthetic cloud field, at various solar zenith angles, illustrate the potential of the method to gain insights into the nature of 3D radiative effects for complicated atmosphere‐surface conditions using information specially derived from the MC simulation. For simulations performed with a single surface albedo it is found that as surface albedo increases, 3D radiative effects increase, too, with maxima occurring at middling to large values, and then decrease. By utilizing the derived coefficients that describe A(x) $\mathcal{A}(x)$ it was established that these 3D effects stem from differences in fractions of radiation entrapped at successive orders of internal multiple reflections for 1D and 3D transfer.

2 of 17 variance reduction techniques can be found in the literature (e.g., Buras and Mayer, 2011;Iwabuchi, 2006;Kutz et al., 2017;Novák et al., 2014;Villefranque et al., 2019;Wald et al., 2014). Since a sampling algorithm will generally behave better if it relies on some knowledge of the processes involved, measuring the efficiency of MC algorithms is also a way to verify that the processes that matter for the quantity of interest (e.g., flux, heating rate, or radiance) are understood well.
Another opportunity offered by MC methods to better understand radiative processes is through examination of simulated paths. In traditional simulations, only the mean and variance of a quantity of interest are output. Although these already contain much information, much more information is contained in the geometry of photon-paths. For example, Várnai and Marshak (2003) analyzed photon-paths to understand how the apparent brightness of a cloudy pixel might be influenced by various parts of a cloud field; Barker et al. (2015) broke estimated radiances into contributions from distinct regions of the domain and various atmospheric constituents in order to aid radiative closure experiments for EarthCARE (Illingworth et al., 2015); Hogan et al. (2019) estimated mean horizontal distance traveled by photons after being scattered to better characterize the entrapment process and evaluate its parametrization. This paper investigates a family of MC methods, referred to as Functionalized Monte Carlo (FMC) methods, and establishes a framework to extract, synthesize, and learn from the large amount of information contained in tracked photon-paths. While standard MC methods output scalar values of variables given specific sets of input parameters (e.g., atmospheric extinction and surface albedos), FMC methods output functional forms of quantities that can be evaluated a-posteriori for input values not used in the simulation.
The idea of FMC methods was branded originally, by William Dunn, as "Inverse MC" because it was used to accelerate the resolution of inverse problems in remote sensing (Dunn, 1981) and medical imaging (Dunn, 1983;Floyd et al., 1986). Recently, Galtier et al. (2017) revisited the FMC methods under the name "Symbolic MC" and showed that they could be extended to output functions of absorption and scattering coefficients for heterogeneous media by combining them with null-collision algorithms (El Hafi et al., 2020;Galtier et al., 2013). Building on these advances, Penazzi et al. (2019) developed an FMC algorithm to solve a coupled conduction, convection, and thermal radiation model in arbitrarily complex geometries and output temperature as a function of convective exchange coefficient. Parallel developments were made by Roger et al. (2019) to extend the method to parameters that were not initially accessible, such as phase function asymmetry parameter, cloud top height, solar angles and viewing geometry.
In the following, the FMC method is used to provide estimates of radiative effects of clouds as a function of underlying surface albedo with a focus on 3D radiative effects of clouds. Studies that have looked into 3D radiative effects of clouds have often focused on dependencies of 3D effects on solar zenith angle and cloud properties such as cloud cover or cloud geometry. For instance, Gristey et al. (2020) investigated the relationship between cloud shape and flux density distribution for highly resolved cloud fields and a number of parameters that capture important properties of the fields. Yet, Hogan et al. (2019) showed that an important mechanism governing cloud 3D radiative effects, in the shortwave, is entrapment of radiation, where a fraction of upwelling radiation is "trapped" by clouds and redirected back to the surface. This fraction is larger when horizontal transport of light is accounted for (i.e., solving for 3D RT rather than the Independent Column Approximation (ICA or 1D)). This is because in 3D, clouds interact with light reflected by the surface in nearby clear-sky areas, but in 1D this light escapes relatively easily to space. While it is known that larger surface albedo directly increases the amount of upwelling radiation that can be entrapped, to our knowledge it has not been explicitly characterized thus far.
In Section 2 the theoretical framework is presented and discussed using single-layer overcast homogeneous clouds. In Section 3, it is used to analyze the dependence of 3D radiative effects of clouds on surface albedo. In Section 4 the scope of the methods, as well as some limitations and prospects, are discussed.

General Principle
In standard MC methods, processes are specified by a set of equations and parameters , and a set of quantities ( ) are estimated. In FMC methods, a standard MC model is run once, with the aim to generate a functional form of A, denoted as ( ) , thereby enabling rapid and accurate, post-simulation, estimation of ( ) for a (wide) range VILLEFRANQUE ET AL.

10.1029/2023MS003674
3 of 17 of . FMC methods might be sorted into two categories: linear and non-linear. In a linear FMC,  is linear with respect to . This is the case whenever the parameter of interest is a source in Green's sense, that is: • a boundary condition, for example, is the solar constant, A is the reflected flux at Top of Atmosphere (TOA).
A single simulation is performed to estimate the domain reflectivity R and ( ) = ; • a volumetric source, for example, is the Planck emission of any layer k, A is the net rate of radiative exchange between two layers i and j. A single simulation is performed to estimate an opticogeometric viewing factor ξ between the two layers, and ( ) = × ( − ) (Eymet et al., 2009); • or an initial condition in transient problems.
An example of a non-linear FMC algorithm is one that simulates high-resolution spectral radiances using a single sample of photon paths. Such algorithms have been implemented before, although not under the name "FMC," for instance for analyzing radiances in the O 2 A-band by Barker (2008, Figure B1) or by Emde et al. (2011) in the MC code MYSTIC (Mayer, 2009). The standard MC approach consists of a single run at wavelength λ using absorption coefficients k a,λ , resulting in, for example, spectral radiances R λ = R(λ; k a,λ ). The FMC approach would, again, perform a single standard MC simulation, now assuming k a = 0, and save a sufficient amount of photon path information to allow subsequent computation of R for any λ across the narrow band. Assuming that the scattering properties of the atmosphere and surface are constant over the spectral band, and having an estimate of the population of photon paths for a gas-free atmosphere, one can compute transmittances along these sampled paths as exp ( − ∫ ( ) ) , which when summed over the path sample leads to quick estimation of high-resolution R λ .
This FMC is facilitated by the fact that atmospheric volume-scattering properties were assumed to be constant in the spectral band, hence the geometry of the paths (number and spatial distribution of successive scattering and surface-reflection events) do not depend on the parameter. Should that no longer be the case, a different FMC could still be designed to account for variations in paths geometry, by further scaling the sampled-path transmissivities, as in Emde et al. (2011). Variation of path geometry with the parameter of interest is one of the main challenges of non-linear FMC methods: if the paths that are sampled during the single simulation are not representative of the paths that would be sampled using other parameter values, then the variance associated with the estimates might be very large. To design an efficient FMC algorithm, one first needs to design a standard MC algorithm in which the geometry of the sampled paths is as independent as possible of the parameter of interest. Once this is achieved, the objective is to factor path information so that: (a) sufficient information is condensed into an amount of data that is manageable; and (b) evaluating a quantity for a given is fast.

A Simple FMC Algorithm: Flux as a Function of Surface Albedo
This section develops an FMC algorithm that gets employed later to illustrate estimation of radiative quantities for cloud fields above reflecting surfaces. More specifically, this FMC uses a standard 3D MC RT model to compute solar cloud radiative effects as a function of surface albedo for a planar, uniform, Lambertian surface. For clarity, it is assumed that underlying surface albedos are spectrally and spatially constant. As discussed in Section 4, however, this is not a fundamental assumption for the method.
In a standard MC simulation run with surface albedo (see Appendix A for a brief description of such a simulation), the flux |̂ is estimated as where |̂ is the weight of the i t h photon, and N is the number of injected photons.
It can be shown that the functional form ( ) of |̂ is a polynomial function of (see for instance the demonstration in Appendix B), where |̂ is the flux after n reflections by a surface of albedo , as illustrated in Figure 1. Note that the leading coefficient, 0|̂ , is atmosphere-only albedo to incident solar spectral irradiance (Barker & Davies, 1992) and is 4 of 17 therefore independent of . In the FMC framework, interpreting Equation 2 as the sum of fluxes reflected n-times by the surface, it can be written as where M is the maximum number of surface reflections undergone by any photon path during the simulation, and the coefficients |̂ are produced by the single FMC simulation by sorting the photon paths by n and letting where N n is the number of photon paths that have undergone n surface reflections. If N n = 0, then so does |̂ , implying that the nth order term does not contribute to |̂ .

Substituting Equation 4 in Equation 3, one can see that
|̂ is a function that scales saved values |̂ , such that where n(i) is number of reflections experienced by the ith path. In this particular "simple" case, all paths that encountered the surface n times receive the same scaling, which results in the condensed expression of Equation 3. The utility of expressing ( ) this way is that individual path weights need not be output. Rather, they can be accumulated during the simulation into counters that correspond to n, thus reducing significantly the amount of data needed to be output by the FMC.
The reason why this FMC is described as being "simple" is because of the level of idealized surface conditions. For a surface with non-uniform albedo, location on the surface of each path encounter would have to be logged. Likewise, if surface albedo depends on incident zenith or azimuth angles, they would have to be logged, too. If, however, a surface is uniform but non-Lambertian with albedo that is independent of incident angles, the FMC developed here applies.
The formal procedure to derive an FMC algorithm from a standard MC algorithm using importance sampling is detailed in Appendix B for the case of the parameter being surface albedo. MC uncertainty associated with the FMC's estimate of ( ) is also a function of α, where whose coefficients can be estimated during the same simulation. The second term ( |̂|̂) involves quantities that are already output from the simulation. To estimate the first term however, more information needs to be recorded during the simulation: one quantity per pair of number of reflections has to be estimated as where ( ) |̂ is the weight of the ith realization after j surface reflections. For j > n(i) (the maximum number of surface reflections of path i), ( ) |̂= 0 . Note that in some MC algorithms, such as the one used here, ( ) |̂ is independent of j for j ≤ n(i), while in other algorithms the weights are modified along the photon paths (e.g., 10.1029/2023MS003674 5 of 17 attenuated due to absorption) and hence may be different for different reflection orders. Note also that for the reflected flux at TOA, ( ) |̂= 0 for all j ≠ n(i), hence only the terms where j = k are non zero.

Setting the Value of T
he essence of the FMC method consists of estimating coefficients of a functional form in which the chosen model parameter(s) is a variable. When FMC algorithms are derived from a standard MC algorithm by reformulating the problem using importance sampling, it introduces arbitrary pdfs that impact the confidence associated with FMC estimates. Indeed, paths that are sampled using a given parameter value during the simulation are not necessarily representative of the path distributions associated with other values of the parameter. Regarding the FMC developed here, the most comprehensive sampling of the population of photon paths is achieved when the standard MC uses = 1 . Indeed, estimates of ( ≤̂) are at least as accurate as the standard MC's ( ) since the paths that would be sampled in this standard MC are subsets of the paths sampled by the FMC. As decreases, however, so too does the ability to sample longer paths, and thus coefficients for the high-order terms in Equation 3 become increasingly uncertain. This results in diminishing quality of estimates of ( ) .
To illustrate this, polynomials (i.e., Equation 3) for surface (SFC) and TOA fluxes as functions of α were estimated, using various values of , for simple atmospheres consisting of homogeneous overcast clouds. Three cases labeled C thin , C med and C thick were studied, each containing a single cloud layer of 1.25 km height and homogeneously filled with liquid water content values of 10 −5 , 10 −4 , and 10 −3 kg kg −1 respectively, resulting in visible optical depths of approximately 2, 20, and 200. Contributions from partially overlapping portions of distributions of photon paths (reflected n times by the surface) to total fluxes, that is, |̂ , are output from the FMC using , for all n ∈ [0, M]. Recall that M, the largest number of reflections undergone by a single path sampled during the simulation, is not chosen a-priori but depends on the atmospheric properties, and the total number of sampled paths N. As any rare event, the probability of sampling at least one path with very large M tends to 1 as N → ∞ (in a scattering media). Figure 2 shows the resulting polynomials along with their uncertainties, as well as polynomials that were truncated at successive orders, and the coefficients output from the FMC simulations. The right panel illustrates that in the simulations using small values of , paths with large numbers of surface reflections, which contribute to the high order terms of the polynomials, are rarely or never sampled. The middle panels shows that if the dependence of the flux on the albedo is quasi-linear, as is the case for the flux reflected by the thin cloud, then the truncated polynomials converge toward the full solution rapidly as the truncature order increases. This means that the FMC estimate will be accurate on the whole [0, 1] range for any value of , even for small for which the longest paths are seldom sampled. However, if the dependence of the flux to the surface albedo is highly non-linear, as is the case for downwelling flux at the surface under optically thick clouds, using small values of and hence poorly sampling the longest paths will lead to inaccurate estimates of the flux at high albedo values (as demonstrated by the gray shading on the red curve of the bottom left plot of Figure 2).
The left panel of Figure 3 shows the relative uncertainty of the downward surface flux estimate associated with (̂) pairs, for the slab of medium optical depth. The relative uncertainty is |̂( )∕|̂( ) and is only defined for 0 . This shows that setting = 1 ensures a minimum uncertainty for all values of α. This remains true for any atmospheric conditions. In the particular case of this FMC algorithm, the value used during the simulation also affects the computing time, as is shown in the right panel of Figure 3: the simulated paths are longer on aver age for larger albedo values , and therefore take more time to be tracked. In the tested cases, simulation time evolves linearly with albedo, and is never more than twice the computation time corresponding to a black surface, whose absolute value depends on the cloud field optical depth. In the presence of thick clouds, the increase of computing time with respect to an increase of surface albedo is slower than for thin clouds. This results from two opposing effects that partially compensate: most paths will be reflected to space before they encounter the surface, yielding null sensitivity of computing time to the value of surface albedo; however, paths that reach the surface will be reflected many times, thereby increasing computing time. Since a single simulation can be used to evaluate fluxes at any surface albedo, it will probably often be worth running a slightly longer simulation to ensure that the longer paths are sampled and the estimated functional can be used with confidence for any value of albedo. In the remainder of this article, all FMC simulations were performed with = 1 .

Using the FMC Method to Learn About the 3D Radiative Effects of Clouds
This section illustrates the broader potential of the FMC method for more complex cloud-radiation problems by analyzing the 3D radiative effects of clouds and their dependency on surface albedo. FMC simulations are conducted on a single cloud field output from a Large-Eddy Model (LEM), as described in the following subsection.

The Cloud Field
The Large-Eddy Simulation used here is an idealized case of shallow cumulus over land, known as the ARM Cumulus case (Brown et al., 2002). It was performed with the French LEM Meso-NH at 8 m spatial resolution on a 12 × 12 × 4 km domain. The LEM was run with a single moment microphysics scheme and no interactive radiation scheme. The latent and sensible heat fluxes at the surface, as well as temperature tendencies due to radiative cooling and large scale advection, were prescribed as a function of time and height.
The 3D field used in the following corresponds to 12:23 Local Solar Time. Figure 4 illustrates the main cloud properties. Cloud cover is defined as the fraction of columns containing any amount of liquid water and equals 38%; when defined as the fraction of columns with optical depth larger than 1, it equals 26%. Cloud base is at 960 m height and the cloud layer is approximately 830 m deep. Fractional standard deviation of in-cloud liquid water mixing ratio (standard deviation normalized by its mean) is 0.62 on average.

How Do 3D Radiative Effects of Cumulus Clouds Depend on Surface Albedo?
The 3D radiative effects on solar downwelling flux at the surface and upwelling at TOA were examined using the algorithm of Section 2. 3D and 1D simulations were performed for eight solar zenith angles (SZA) from 0 (at zenith) to 75°. From these simulations, the polynomial coefficients for surface and TOA fluxes were output, along with their covariances from which the statistical uncertainties around the mean polynomial functions, themselves polynomial functions of the albedo, were constructed according to Equations 6 and 7. By means of an extremely fast post-treatment (evaluating a polynomial function), estimates of surface and TOA fluxes, as well as their uncertainties, were estimated for 200 values of surface albedo regularly spaced between 0 and 1. The 3D effects are differences between flux estimates from the 3D and 1D simulations. The RT simulations were run with constant effective radius of 10 μm. Assuming that droplets sizes tend to increase going from edges of clouds to their cores, accounting for variations in effective radius would tend to suppress variations in optical depth, which would shift 3D RT results slightly toward 1D results. Although this remains to be verified, it is assumed that the shape of the functional would only marginally be affected, while the amplitude of 3D effects would be slightly smaller. Figure 5 presents the 3D effects at surface and TOA as functions of surface albedo. First, for the surface, over black and dark surfaces (small albedo), 3D effects are positive when the sun is high (small SZA) and negative when it is low (large SZA). This means that at small SZA, 3D effects increase surface downwelling flux, while the opposite is true at large SZA. The effect is opposite at TOA (less reflective clouds at small SZA and more reflective clouds at large SZA). This effect and the underlying mechanisms have been detailed and described many times illustrating the competition between partially compensating effects that depend on SZA. For 3D calcula tions, direct radiation from the sun is intercepted by cloud sides and is more likely reflected to space which makes larger shadows and brighter clouds than in 1D. In 3D the downward flux at surface is affected by more diffuse radiation due to photons from direct beam being reflected downward by cloud sides, downward-traveling photons escaping through cloud sides, and enhanced entrapment, that is, upward-traveling photons being reflected by surrounding clouds back toward the surface (e.g., Hogan et al., 2019;Várnai and Davies, 1999). Second, as surface albedo increases toward medium values, so do 3D effects at the surface, meaning that clouds bring more flux to the surface in 3D than in 1D for all SZAs, while the 3D effects at TOA decrease, meaning that clouds become less and less reflective in 3D compared to 1D. Third, 3D effects at surface reach a maximum (a minimum at TOA) for a value of albedo that depends on SZA and is not the same at the surface and TOA; thereafter they start to decrease at the surface (increase at TOA).
To better understand this behavior, the difference between the 3D and 1D polynomial coefficients obtained with = 1 , which are the contributions of the different orders of reflection into the total 3D effects at α = 1, are presented in Figure 6. The 3D and 1D coefficients are also plotted separately. The contribution of the nth order of reflection to the same quantity at another albedo α is simply A n|1 × α n . The zeroth-order of reflection corresponds to fluxes at α = 0.
As previously stated, the zeroth-order of 3D effects strongly depends on SZA with either positive or negative effects. The sign of the once-reflected 3D effects, however, no longer depends on SZA: cloud fields in 3D systematically bring more flux to the surface and less to TOA than in 1D after one surface reflection. This is a signature of 3D effects on entrapment, as 3D photons that reach the surface via clear-skies can be trapped by neighboring clouds and redirected back to the surface, while for 1D they stand a much better chance of escaping to space. In Figure 5, this translates into the positive (negative) slope seen in the surface (TOA) 3D effects at small albedos. The magnitude of the coefficient (and hence steepness of the line) does, however, depend on SZA, mostly because it is a fraction of the flux that reached the surface at the zeroth-order, which is itself affected by SZA. After a second reflection, 3D effects become negative at the surface and positive at TOA. This is because 10.1029/2023MS003674 9 of 17 3D photons that were trapped by clouds after one reflection managed to escape to space after a second reflection, whereas the 1D photons that were trapped under the clouds are still there.
Looking at the 1D and 3D coefficients independently shows that the 3D coefficients decrease much more rapidly with the number of surface reflections that the 1D coefficients. Let us write the nth order surface flux A n as an atmosphere-reflected fraction r n of the surface-reflected fraction of A n−1 such that when α = 1, This is often used in operational RT models with the assumption that r n = r for all n, as is detailed in Appendix C.
In the bottom panels of Figure 6, a constant r n for all n would yield a straight line because of the log scale of the y axis. The fact that the 3D coefficients are closer to a straight line than the 1D coefficients implies that after each surface reflection, a similar fraction of the 3D fluxes is reflected by the atmosphere back to the surface, whereas this fraction increases with n for 1D fluxes. Figure 7 shows maps of A n coefficients for downwelling flux at the surface for the first 3 orders n and = 1 , as well as the second order structure functions of the A n maps as a function of scale for the first 5 orders (Barker et al., 2017). These coefficients were obtained from additional simulations performed with a slightly different 3D MC model that is optimized for spatially discretized radiative estimates such as flux maps and physically-based rendering (htrdr, Villefranque et al., 2019).
where A is a two dimensional square map of shape (N x , N x ), Δx is the horizontal resolution of the map and k an integer in [1, N x ]. S 2 (L) is large when the amplitude of the field is large. At small L, S 2 measures the average variability inside small regions (e.g., inside a cloud or a clear sky hole). As L increases, so does S 2 if the field is structured into patterns of scale smaller than L. If the field is smooth, then S 2 is independent of L. Figure 7 shows that the fields are characterized by patterns of approximately 1 or 2 km for the zeroth-order in 3D and all orders in 1D. For subsequent orders in 3D, quasi-uniform maps and horizontal curves demonstrate that the radiative field is almost entirely smoothed by horizontal transport after the first surface reflection. The structure of the 1D fields however is largely unchanged after all reflections, confirming that the relatively large coefficients for high reflection orders shown in Figure 6 are mostly due to photons that are trapped within cloudy columns and continue to bring light to the surface reflection after reflection.

Summary and Discussion
Analysis of the optical paths space was performed using a family of MC algorithms that output functionals of chosen variables instead of scalar values. This has allowed us to study the dependency of cloud radiative effects on surface albedo from one unique MC simulation instead of many simulations run with different surface albedo values. It was shown that 3D radiative effects of clouds reach a maximum for medium to large albedo values because of the increasing intensity of photon entrapment, and then decreases as more entrapped 3D radiation manages to find clear sky holes through which to escape to space, somewhat compensating the 1D overestimation of cloud reflectivity. After a few reflections at the surface there are small amounts of entrapped 3D radiation while there remains 1D radiation trapped under clouds, enhancing surface illumination, especially over very bright surfaces. Note that a single cloud field was used in this study, hence the sensitivity to these conclusions to cloud characteristics such as cloud cover, liquid water path, cloud base height or cloud geometrical thickness was not investigated. From the literature and the few sensitivity tests that were performed, it can be assumed that the amplitude of 3D effects increases with cloud geometrical thickness and total cloud cover up to medium values (e.g., Črnivec & Mayer, 2019;Gristey et al., 2020). In overcast cloud fields, 3D effects are generally less intense due to the lack of cloud sides for radiation to interact with, which might modify the shape of the functional in addition to its amplitude.
Beyond these findings, the main message here is that MC methods are much more than a way to compute accurate fluxes. There is a very large potential for learning from the paths sampled during the simulation. The Functionalized MC method provides a theoretical framework in which this can be achieved. Designing a FMC algorithm will often require a shift in perspective on the object of study; for instance here seeing the surface flux as the sum of n-times reflected fluxes. This renewed perspective naturally enhances the power of analysis already provided by the simulations. FMC methods are also a way to synthesize the information contained in the path samples, which might otherwise be very difficult to comprehend or even technically impossible to manipulate.
It was our intention to emphasize this aspect rather than the perhaps more obvious benefits of the FMC methods: the fact that they greatly accelerate MC methods, as only one simulation is needed to recover the amount of information that would normally be output from an arbitrary large number of standard simulations. This property has been useful in particular in the context of parameter identification or inversion (Dunn, 1981;Floyd et al., 1986;Maanane et al., 2020;Sans, Blanco, et al., 2021) and more generally has a great potential for global nonlinear parametric sensitivity analyses. In other words, FMC methods can be thought of as a physically-informed machine learning technique.
Surface albedo is a parameter that is very well suited to FMC methods because the level of data compression can be high, and because an optimal value can be used for importance sampling ensuring that small uncertainties can always be achieved for the whole parameter definition domain. A small amount of information was saved from the RT simulations due to the assumption of α being spectrally and spatially constant. Should this assumption be relaxed, additional information would be required to reconstruct the functional. The most straightforward way would be to store the location and wavelength of every surface reflection. Then, given a field of spectrally and spatially varying surface albedo, sampled path weights would be corrected accordingly, preserving correlations between spatial and spectral variations of atmosphere and surface.
The fact that = 1 is optimal is due to the fact that paths that would be sampled with smaller are actually included in the paths sampled with = 1 ; these longer paths only need to be truncated after a few reflections to recover the shorter paths without information loss. This is not always the case and this is one of the limitations of importance sampling based FMC methods. The recent work of Roger et al. (2019) and Penazzi et al. (2019) has extended FMC methods beyond use of importance sampling, to cases where either the parameters of interest do not appear explicitly in the RT equation or the variance cannot be controlled; at the moment it remains limited to a few examples such as solar zenith angle and phase function asymmetry parameter, but the potential is huge. Another current limitation of the method (which is not related to importance sampling) is the difficulty to handle large numbers of parameters or more complex parameters (such as the whole field of extinction or effective radius), because of the increased amount of data needed to be stored, and the complexity of the resulting functional, implying that a significant amount of time might become necessary to evaluate the functional for a given parameter set.
Finally, it seems that there is still much to discover on the properties of FMC methods; further investigations will help deepen our understanding of these methods and of their significance, exploring their potential and their links with other statistical methods, and extending them to yet inaccessible configurations. With all this, and although still a long-term perspective, it does no longer seem inconceivable that MC-based methods could be used to solve for 3D RT in atmospheric models.
To estimate A with a Monte Carlo (MC) algorithm, first a wavelength λ and an entry point (x, y) are sampled and then a path is tracked to estimate F. A path is defined as a list of vertices, each consisting of a location x = (x, y, z) and direction ω. The first vertex of any path that contributes to F is x 0 , ω 0 , and its last vertex is located at z = z TOA and in any upwelling direction ω ↑ . Between these first and last, any number of vertices can exist, in particular any number of vertices located at the surface and hence any number of surface reflections. A path can be defined as a random variable Ψ on a path space Ω Ψ with probability density function p Ψ (ψ). Then if the contribution of a path Ψ to F is w(Ψ), This is a conceptual, condensed form representing multiple scattering and absorption in a heterogeneous medium and surface interactions. It could be expanded into a more classical notation where the successive events would appear explicitly; this was not done here for the sake of brevity. Because surface reflections are the focus of this article, let us write F in a form that keeps the scattering sub-paths in a condensed format but lets surface interactions appear explicitly, as where Γ is a "multiple-scattering sub-path" random variable characterized by the path space Ω Γ and its probability density function p Γ (γ|x, ω). These paths are similar to Ψ paths with the constraints that they start and end either at the surface or at TOA, with all other vertices (if any) located inside the domain. In other words, they are the (multiply scattered) sub-paths that connect surface to surface, surface to TOA or TOA to TOA. If a sub-path starts at x, ω then its last vertex is noted x′, ω′. H(x′ ∈ toa) = 1 if γ reaches TOA and 0 otherwise (which means that γ has ended at the surface). F ⊙ is incoming solar radiation at TOA. ≡ ( ′ , ′ ) is surface albedo seen by 13 of 17 sub-path i, which in all generality may depend on wavelength, location and direction of the path when it reaches the surface. p Ω (ω j |ω i ) is the normalized surface reflectivity function that gives the probability of the path being reflected into ω j from incoming direction ω i . Upon reaching the ground, the photon gets reflected with probability α i while upon reaching TOA, the path ends and its contribution is added to the flux counter. Other paths (absorbed by the medium or surface) do not contribute to reflected flux at TOA hence their contributions are zero and they need not appear explicitly in Equation B3. Other formulations can be derived for other quantities such as surface flux or layer absorption.
Equation B3 can be factorized into where F n is the contribution to the total reflected flux by paths that have survived n ≥ 0 surface reflections. Therefore, Then the parameters α i can be isolated from the rest of the expression using importance sampling (in which all the following expressions are for n > 0), yielding where are positive albedo values chosen arbitrarily. If and α i are independent of the incoming direction and uniform in space and spectrum, that is, independent of the path properties (which is not necessary but simplifies the expression here), then Recognizing that the term following (̂) is similar to F n from Equation B6 but with instead of α, and so The terms |̂ are of the same nature as the quantity that would be estimated by the standard MC algorithm of Equation B3, but they result from partial integration over only a subsample of the path space and not the whole path space, and correspond to an albedo of instead of α.
The more general steps to derive a FMC algorithm from a standard MC algorithm using importance sampling are summarized in Table B1.
14 of 17 Step   Table B1 Steps to Derive a Functionalized Monte Carlo Algorithm From a Monte Carlo Algorithm