A Method for Estimating the Cloud Adjacency Effect on the Ground Surface Reflectance Reconstruction from Passive Satellite Observations through Gaps in Cloud Fields

Abstract

A method for estimating the cloud adjacency effect on the reflectance of ground surface areas reconstructed from passive satellite observations through gaps in cloud fields is proposed. The method allows one to estimate gaps of cloud fields in which the cloud adjacency effect can be considered small (the increment of the reflectance $\Delta r_{surf}\le$ 0.005). The algorithm is based on statistical simulation by the Monte Carlo method of radiation transfer in stochastic broken cloudiness with a deterministic cylindrical gap. An interpolation formula is obtained for the radius of the cloud adjacency effect that can be used for the reconstruction the ground surface reflectance in real time without calculations by the Monte Carlo method.

1. Introduction

The reflectance of the ground surface is one of the main surface parameters reconstructed from satellite data. It is widely used for the monitoring of the vegetative cover of the ground surface, the state of water reservoirs, and other objects. A number of algorithms for reconstructing the ground surface reflectance from satellite data have been developed for cloudless situations [1,2,3,4,5]. However, if the ground surface is covered by low-lying clouds, the reflectance cannot be reconstructed from satellite measurements. In addition, images of cloudless areas are also significantly affected by the cloud adjacency effect. In this regard, the problem arises of estimating the radius of the cloud adjacency effect for which the atmosphere can be considered cloudless, and the clear sky algorithms [1,2,3,4,5] can be used for reconstruction of the ground surface reflectance.

This problem was studied in a number of works (e.g., [6,7,8,9,10,11,12]). Cahalan et al. [6] estimated the influence of a single cloud shaped as a parallelepiped on the results of reconstruction of the reflectance of a cloudless area on the ground surface. It was shown that for observations in the nadir in the Landsat B1 channel for the aerosol optical depth (AOD) of the cloudless atmosphere equal to 0.2 and the optical thickness of cloudiness equal to 20, this effect increased the reflectance $r_{surf}$ by 0.015 on the sunward side of the principal plane and decreased $r_{surf}$ by 0.03 on the antisun side of the principal plane. Away from the cloud, its influence on $\Delta r_{surf}$ decreased and reached positive asymptotic values less than 0.005 at a distance of 2–3 km. In the B7 Landsat channel for the AOD of the cloudless atmosphere equal to 0 and the cloud optical thickness equal to 20, this influence increased $r_{surf}$ by 0.004 on the sunward side of the principal plane and reduced $r_{surf}$ by 0.048 on the antisun side of the principal plane. Away from the cloud, its influence on $\Delta r_{surf}$ decreased and vanished at a distance of 3 km. The main disadvantage of this approach is that the estimate did not take into account the interaction of clouds in cloud fields.

Nikolaeva et al. [7] considered continuous half-plane cloudiness. It was shown that in some situations, the effect of cloudiness on the accuracy of reconstruction of the ground surface reflectance could not be neglected at distances from neighboring clouds no less than 25 km. However, the cloud field is broken in most real situations, and the radius of the cloud adjacency effect can differ significantly.

Wen et al. [8] studied two test scenes with pixels screened by broken cloudiness. For a deterministic cloudiness model, the dependence of the reflectance on the radius of the cloud adjacency effect was calculated by the Monte Carlo method for particular arrangements of clouds and optical and geometrical conditions of observations. It was found that for scene 1, at a distance of 4 km from cloudiness, the reflectance increased by around 0.01 at a wavelength of 0.47 $\mathsf{\mu }$m and by around 0.003 at a wavelength of 0.66 $\mathsf{\mu }$m. For scene 2 at a distance of 8 km, it increased by around 0.006 at a wavelength of 0.47 $\mathsf{\mu }$m and by around 0.004 at a wavelength of 0.66 $\mathsf{\mu }$m. Statistical simulation of the radiation propagation in complex three-dimensional media by the Monte Carlo method (analogous to [8]) for real-time solution of the inverse problems of passive ground surface sensing from space is inexpedient because it is very labor-intensive and the results obtained cannot characterize all possible situations due to their peculiarities.

Va’rnai et al. [9] used statistical analysis of images of the northern region of the Atlantic Ocean obtained with the MODIS spectroradiometer in 2000–2007 to study the adjacency effect. Average reflectance was obtained in five MODIS channels as functions of the adjacency effect radius and their increments $\Delta r_{surf}$ were estimated. On average, the adjacency effect was observed for radii up to 15 km. The effect of cloudiness on the ground surface reflectance reconstructed near sunlit and shadowy sides of clouds differed significantly for radii up to 4 km. For larger distances, the reflectance weakly depended on the cloud arrangement. This approach has the following restrictions: (1) large variance of the results obtained caused by wide variability of the optical and geometrical conditions; (2) results were obtained only for the ground surface covered with water, and their extrapolation to the land surface was complicated by high horizontal surface inhomogeneity.

Marshak et al. [10] considered a uniform model stochastic field of parallelepiped clouds with identical sizes for the nonreflecting ground surface. The cloud adjacency effect averaged over the cloud field was estimated. It was reported that for the nonreflecting surface, the result obtained weakly depended on the cloud size but strongly depended on the cloud cover index ${\delta }_{cl}$. The effect of cloudiness on images of cloudless areas was estimated as a function of the optical cloud thickness and the cloud cover index. It was demonstrated that the reflectance can increase by $\Delta r_{surf}\sim$ 0.1 for high cloud cover indices (${\delta }_{cl}\sim$ 1) and optical cloud thicknesses (${\tau }_{cl}\sim$ 50).

Marshak et al. [11] showed that for weakly reflecting surfaces (with the reflectance close to 0), the adjacency effect decreases linearly with increasing wavelength in the visible range; therefore, in these situations, it was sufficient to consider the smallest radiation wavelength. The algorithm was constructed as follows:

• The reflectance for a cloudless area $r_{surf,clear}$(${\lambda }_1$) was determined at the smallest wavelength (${\lambda }_1=0.466$ $\mathsf{\mu }$m) under the assumption that only molecular scattering occurs beyond the cloudiness.
• The parameters a and b that relate the reflectance of areas strongly affected by cloudiness at the wavelengths ${\lambda }_1$ = 466 nm and ${\lambda }_2$ = 855 nm ($r_{surf,cloud}\left({\lambda }_2\right)=a{r}_{surf,cloud}\left({\lambda }_1\right)+b$) were determined from measurements in the 20 × 20 pixel window.
• The reflectance of the area shadowed by cloudiness at the wavelength ${\lambda }_2$ was corrected for the average adjacency effect at the wavelength ${\lambda }_2$ obtained by multiplication of the average adjacency effect at the wavelength ${\lambda }_1$ to the constant a, i.e.,

  $$r_{surf,clear}\left({\lambda }_2\right)=r_{surf,cloud}\left({\lambda }_2\right)+a\Delta r\left({\lambda }_1\right),$$

  (1)

  where $r_{surf,cloud}\left({\lambda }_2\right)$ is the ground surface reflectance at the wavelength ${\lambda }_2$ changed in the presence of cloudiness, and $\Delta r\left({\lambda }_1\right)$ is the increment of the reflectance in the presence of cloudiness at the wavelength ${\lambda }_1$ averaged over the cloudy field.

In the above paper [11], the possibility of application of this approach was shown for a wide range of model optical and geometrical conditions. However, the approaches developed in papers [10,11] can be inapplicable if the surface is not weakly reflecting.

In our work [12], a model of continuous cloudiness with a deterministic cylindrical gap and the ground surface reflectance changing from 0 to 1 was considered. Results of our calculations showed that when the gap radius R increased from 4 to 15 km, the recorded radiation intensity changed by less than 10% of its value under conditions of clear sky depending on the cloudiness. In the present work, the problem formulation [12] discussed above is used, but under conditions of broken cloudiness.

The key to solving this problem is the choice of a cloud field model. For these purposes, continuous [7,12], deterministic [8,11], Gaussian [13,14,15], and Poisson models of cloudiness can potentially be used with parallelepiped [10,13,16] or paraboloid clouds [13,14]. Below, we use the Poisson models of broken cloudiness with paraboloid clouds. The reason is that these models fit well for statistical simulation of radiative transfer in broken cumulus cloudiness and do not require large computational time for the generation of individual realizations of cloud fields. In the present article, we deliberately do not consider the deterministic models of cloud fields because, in this case, the results obtained will be rigidly adhered to a specific realization of the cloud field.

Radiation transfer in the broken cloud field can potentially be studied by the following methods. In papers by Wen et al. [8] and Varnai et al. [9], the Monte Carlo method was used for statistical simulation of radiation transfer in a three-dimensional inhomogeneous deterministic medium. Some authors [10,14,15,16,17] have used the Monte Carlo method with averaging of the results over an ensemble of realizations of a cloudy field for statistical simulation of radiation transfer in broken cloudiness. Zuev and Titov [13] and Titov et al. [18] used the closed-form equation method, according to which the problem was solved for the effective horizontally homogeneous medium such that the received radiation intensity coincided with the radiation intensity averaged over an ensemble of realizations of the cloud field. This method is most effective from the viewpoint of the computing time. However, as shown in papers [13,19], the closed-form equation method has restrictions on the solar zenith angles and the satellite zenith angle. Below, we use the Monte Carlo method with averaging over the ensemble of cloud field realizations.

Our approach presented below is free of all the above-enumerated disadvantages of the alternative approaches.

2. Problem Formulation and Solution Method

The problem is formulated as follows (Figure 1).

Let there be the flat atmosphere—the ground surface system. The parallel solar radiation flux is incident on the upper boundary of the atmosphere at the angles ${\theta }_{sun}$. Let the statistically homogeneous Poisson cloud field consist of paraboloids. The algorithm of constructing the Poisson fields of paraboloids described in [13] was used for simulation. The number of clouds was random and obeyed the Poisson distribution with mathematical expectation presented in [13] (pp. 119–125). The positions of cloud centers were random and uniformly distributed in the horizontal plane. The cloud shapes were similar, and their sizes were random and obeyed the exponential distribution [13] (pp. 119–125) with the mathematical expectation equal to the average horizontal cloud size $\overline{L}$. The deterministic cylindrical gap of radius R was simulated within the cloud layer. The lower boundary of cloudiness was fixed and equal to $h_{min}$. The thickness of a separate cloud was random, but the average thickness was equal to $\Delta h$.

Let the cloud cover index be equal to ${\delta }_{cl}$. The aerosol optical characteristics were homogeneous within the cloud medium and were determined by extinction and scattering coefficients and a scattering phase function. The optical model of cloudiness was constructed based on the OPAC cumulus cloud model [20]. The optical characteristics of the cloudless atmosphere were described by the LOWTRAN-7 model [21]. The homogeneous Lambert ground surface with the reflectance $\overline{r_{surf}}$ was considered. The satellite sensing system was placed at the altitude $h_d$. Its optical axis was directed toward the point on the ground surface located in the center of the projection of the gap in the cloud field onto it. Radiation at the wavelength $\lambda$ was received at this point. The satellite zenith angle was ${\theta }_d$, and the azimuthal angle between the directions toward the receiver and the Sun from the observation point on the ground surface was equal to $\phi$.

It was necessary to estimate the gap radius $R_{*}$ for which the neglect of the cloud adjacency effect will introduce the error in the reconstructed reflectance $\Delta r_{surf}$ less than 0.005 and will allow the mask of pixels to be constructed for which the reconstruction error $\Delta r_{surf}$ will exceed 0.005.

The problem is solved in several steps:

(1)

From the MODIS data, the cloudiness mask (MOD06L2 data), the AOD of cloudless areas (MOD04L2), the upper cloud boundaries (MOD06L2), the AOD of cloudiness (MOD06L2), the solar zenith angles ${\theta }_{sun}$, and the satellite zenith angle ${\theta }_d$ and the azimuth angle $\phi$ of the receiving system (MOD03L2) are determined.

(2)

The AOD value averaged over clear sky pixels is determined.

(3)

Among the LOWTRAN-7 cloudless models, the model is selected with the AOD value closest to that of the MODIS channel.

(4)

From the MOD09L2 data, the average reflectance $\overline{r_{surf}}$ of cloudless areas is determined.

(5)

From the AOD of cloudiness and the altitude of the upper cloud boundary $h_{max}$, the cloud extinction coefficient is determined. The cloud medium is considered homogeneous. The quantum survival probability and the normalized cloud scattering phase function are selected from the LOWTRAN-7 models.

(6)

The molecular scattering coefficient is selected from the mid-latitude summer LOWTRAN-7 models.

(7)

The radiation intensities $I_{sum,cloud}\left(R\right)$ received by the satellite system and averaged over an ensemble of realizations of the cloud field for the examined optical and geometrical conditions, the preset average reflectance $\overline{r_{surf}}$, are calculated by the Monte Carlo method depending on the gap radii R.

(8)

The cloud adjacency effect on the reconstructed ground surface reflectance is estimated based on the expression for the total radiation intensity received by the satellite system in the cloudless case for the homogeneous ground surface (in the independent pixel approximation):

$$I_{sum,clear}=I_{sun,clear}+\frac{r_{surf}{E}_{0,clear}}{1-r_{surf}{\gamma }_{1,clear}}{\tilde{I}}_{surf,clear},$$

(2)

where $I_{sun,clear}$ is the radiation intensity scattered in the atmosphere and not interacted with the ground surface for the cloudless atmosphere; $E_{0,clear}$ is the sun irradiance of the ground surface disregarding reflected radiation for the cloudless atmosphere; ${\gamma }_{1,clear}$ is the atmospheric albedo (contribution of singly reflected radiation to the ground surface illumination in the cloudless atmosphere); ${\tilde{I}}_{surf,clear}$ is the intensity of radiation reflected from the ground surface and recorded with the receiving system for the unit ground surface luminosity of the cloudless atmosphere.

From Formula (2), we obtain the expression for the reflectance $r_{surf}$:

$$r_{surf}=\frac{Q⁄E_{0,clear}}{1+{\gamma }_{1,clear}Q⁄E_{0,clear}}$$

(3)

where

$$Q=\frac{I_{sum,clear}-I_{sun,clear}}{{\tilde{I}}_{surf,clear}}.$$

(4)

If we neglect the cloud adjacency effect at the point on the ground surface corresponding to the projection of the center in the gap of the cloudy field, we obtain the approximate value of the reflectance ${\tilde{r}}_{surf}$:

$${\tilde{r}}_{surf}\left(R\right)=\frac{\tilde{Q}⁄E_{0,clear}}{1+{\gamma }_{1,clear},\tilde{Q}⁄E_{0,clear}}$$

(5)

where

$$\tilde{Q}=\frac{I_{sum,cloud}\left(R\right)-I_{sun,clear}}{{\tilde{I}}_{surf,clear}},$$

(6)

$I_{sum,cloud}\left(R\right)$ is the total radiation received for satellite observations of the point on the ground surface located in the center of the deterministic gap with radius R.

Then, the gap radius $R_{*}$ for which the neglect of the cloud adjacency effect introduces the reflectance reconstruction error $\Delta r_{surf}$ less than 0.005 is defined as the least radius $R=R_{*}$ for which the condition

$$\Delta r_{surf}=\mid r_{surf}-{\tilde{r}}_{surf}\left(R_{*}\right)\mid \le \delta =0.005$$

(7)

is satisfied.

(9)

Proceeding from the radii $R_{*}$, the mask of pixels is constructed for which the cloud adjacency effect is significant.

The radiation intensity $I_{sum,cloud}\left(R\right)$ received by the satellite system and averaged over an ensemble of realizations of the cloud field was simulated by the backward Monte Carlo method [22]. The trajectories were subdivided into P packages, and each package comprised N trajectories. The mean free path length beyond the cloud layer was simulated by the standard algorithm [22] (p. 11). For the photon trajectories passing through the cloud layer $[h_{min},h_{min}+\Delta h]$, the maximum cross-section method [22] (p. 12) is widely used to simulate the mean free path length. However, our estimations have shown that the standard algorithm [22] (p. 11) adapted by us for computations in a three-dimensional inhomogeneous medium is more economic from the viewpoint of the computational time. The scattering and absorption of photons in the medium was simulated by standard algorithms [22] (p. 10). When a photon collided with the ground surface, its reflection was simulated, and its weight decreased by the probability of radiation absorption by the surface. At each collision point of the photon trajectory in the medium and collision with the ground surface, the local estimate of radiation coming at the upper boundary of the atmosphere in the direction toward the Sun was calculated from the formulas:

For the atmosphere

$$I_{pnk}=S_{\lambda }{q}_{pnk}\frac{{\sigma }_{sa,pnk}{g}_{a,pnk}\left({\mu }_{pnk}\right)+{\sigma }_{sm,pnk}{g}_{m,pnk}\left({\mu }_{pnk}\right)}{{\sigma }_{t,pnk}}exp(-{\tau }_{pnk}),$$

(8)

and for the ground surface

$$I_{pnk}=S_{\lambda }{q}_{pnk}\frac{r_{surf}}{\pi }exp(-{\tau }_{pnk}),$$

(9)

where $I_{pnk}$ is the estimated intensity of radiation coming at the upper boundary of the atmosphere from the kth collision point of the nth trajectory from the pth package; $q_{pnk}$ is the energy (weight) of the photon at the collision point; $S_{\lambda }$ is the solar constant, ${\sigma }_{sa,pnk}$ is the aerosol scattering coefficient, ${\sigma }_{sm,pnk}$ is the molecular scattering coefficient, and ${\sigma }_{t,pnk}$ is the extinction coefficient at the collision point, $g_{a,pnk}\left({\mu }_{pnk}\right)$ is the phase function of aerosol scattering at the angle ${\mu }_{pnk}$, $g_{m,pnk}\left({\mu }_{pnk}\right)$ is the molecular scattering phase function, and ${\tau }_{pnk}$ is the optical thickness of the layer from the collision point to the upper boundary of the atmosphere in the direction opposite to the direction of radiation incidence at the upper boundary of the atmosphere.

If the photon interacted with the medium within the cloud layer, the cloud optical aerosol characteristics were considered in Formulas (8) and (9), and if beyond the cloud layer, the aerosol characteristics of the cloudless layer were considered. It is difficult to calculate ${\tau }_{pnk}$ from Formulas (8) and (9), because the collision point can lie below, within, or above the cloud layer; the ray trajectory from the collision point toward the Sun can intersect the deterministic gap of the cloud field; and it also can intersect cloud boundaries or clouds can partially overlap each other.

Then, the total radiation intensity is calculated from the formula:

$$I=\frac{1}{PN}\sum _{p=1}^P\sum _{n=1}^N\sum _{k=1}^{K_{pn}}{I}_{pnk},$$

(10)

where P is the number of trajectory packages, N is the number of photon trajectories in one package ($P\times N$ trajectories altogether), $K_{pn}$ is the number of collisions for the nth trajectory from the pth package. The error of estimate (10) is calculated from the formula:

$$\delta I=\frac{\sqrt{\frac{1}{PN}{\sum }_{p=1}^P{\sum }_{n=1}^N{\sum }_{k=1}^{K_{pn}}{\left(I_{pnk}\right)}^2-{\left(\frac{1}{PN}{\sum }_{p=1}^P{\sum }_{n=1}^N{\sum }_{k=1}^{K_{pn}}{I}_{pnk}\right)}^2}}{\frac{1}{PN}{\sum }_{p=1}^P{\sum }_{n=1}^N{\sum }_{k=1}^{K_{pn}}{I}_{pnk}}.$$

(11)

3. Testing of the Algorithm

The results of calculations of $I_{sum,cloud}$ (0.5) for stochastic paraboloid clouds were compared with those for equidistant deterministic paraboloids in the situation with $\lambda =0.55$ $\mathsf{\mu }$m, R = 0.5 km, L = 1 km, ${\theta }_{sun}=0^{\circ }$, ${\theta }_d$ = 0 and 45${}^{\circ }$ ${\tau }_{cl}$ = 20, ${\tau }_0$ = 0.285, ${\delta }_{cl}$ = 0.5, $h_{min}$ = 0 km, and $\Delta h$ = 4 km (Figure 2). The results of test calculations are shown in Figure 3. Calculations were carried out for P = 30 and N = 5000. We obtained $\delta I\le$ 0.004 for the equidistant paraboloids and $\delta I\le$ 0.034 for averaging over the ensemble of random paraboloid realizations.

The results of our comparison showed that $I_{sum}$, cloud for equidistant clouds was less than $I_{sum}$, cloud averaged over an ensemble of cloud field realizations (on average, by 9.7% in Figure 3a and by 12.8% in Figure 3b). Nevertheless, the results were close to each other, which allows us to state that the developed program with averaging over an ensemble of realizations of the cloud field works correctly.

4. Interpolation Formula for Estimating ${\mathit{R}}_{*}$

Using the developed program of the Monte Carlo method, we calculated $R_{*}$ values for five MODIS channels and a wide range of optical and geometrical conditions. Calculations were carried out for the following conditions:

• 5 MODIS channels ($\lambda$ = 0.41, 0.47, 0.55, 0.68, and 0.86 $\mathsf{\mu }$m);
• Cloud layer thicknesses $\Delta h$ = 0.5, 1.5, 4, and 6 km;
• Cloud cover indices ${\delta }_{cl}$ = 0.05, 0.15, 0.3, and 0.5;
• Solar zenith angles ${\theta }_{sun}$ = 0, 30, 45, and 60${}^{\circ }$;
• Zenith angles of the receiving system ${\theta }_d$ = 0, 30, 45, and 60${}^{\circ }$;
• Azimuthal angle $\phi =0^{\circ }$;
• Ground surface reflectances $r_{surf}$ = 0, 0.1, 0.3, 0.5, and 1;
• Midlatitude summer model;
• Aerosol optical depth of the cloudless atmosphere ${\tau }_{0.55}$ = 0, 0.09, 0.3, and 0.89;
• Cloud extinction coefficients ${\sigma }_{cl}$ = 10, 20, 30, and 40 km${}^{-1}$.

Calculations were carried out only for one azimuthal angle $\phi$ because our estimations showed that the cloud adjacency effect is maximal at $\phi$ = 0${}^{\circ }$. Calculations were carried out for 102,400 computational grid nodes.

Based on the results of calculations, the program was developed that allowed the radius of the cloud adjacency effect to be estimated, obviating the necessity of simulation of radiation transfer in an atmospheric layer. The program is given in [23]. In the program, generalization of bilinear interpolation over the seven-dimensional grid of values is used:

$$\begin{array}{c}{R}_{*}=C_0+C_1\Delta h+C_2{\delta }_{cl}+C_3(1-{\mu }_{sun})+C_4(1-{\mu }_d)+C_5{r}_{surf}\hfill \\ +C_{6}exp(-\tau )+C_7{\sigma }_{cl}+C_{1,2}\Delta h{\delta }_{cl}+\cdots +C_{6,7}exp(-\tau ){\sigma }_{cl}\\ +C_{1,2,3}\Delta h{\delta }_{cl}(1-{\mu }_{sun})+\cdots +C_{5,6,7}{r}_{surf}exp(-\tau ){\sigma }_{cl}+\cdots \\ \hfill +C_{1,2,3,4,5,6,7}\Delta h{\delta }_{cl}(1-{\mu }_{sun})(1-{\mu }_d)r_{surf}exp(-\tau ){\sigma }_{cl}\end{array}$$

(12)

where $C_0$, $C_1$, ⋯, $C_{1,2,3,4,5,6,7}$ are constants adjusted so that the interpolation formula yielded results coinciding in the computational grid nodes with the calculated ones; $\tau$ is the AOD of the cloudless atmosphere for the corresponding channel.

To estimate the possibilities of the interpolation formula, two series of comparisons were performed for the following optical and geometrical conditions: (1) $\lambda$ = 0.41, 0.55, and 0.86 $\mathsf{\mu }$m, ${\tau }_{0.55}$ = 0.09, ${\delta }_{cl}$ = 0–0.5, $\Delta h$ = 1.5 km, ${\theta }_{sun}={\theta }_d={45}^{\circ }$, $\phi =0^{\circ }$, ${\sigma }_{cl}$ = 20 km${}^{-1}$, and $r_{surf}$ = 0.1; (2) $\lambda$ = 0.41, 0.55, 0.86 $\mathsf{\mu }$m, ${\tau }_{0.55}$ = 0.09, ${\delta }_{cl}$ = 0.3, $\Delta h$ = 1.5 km, ${\theta }_{sun}={\theta }_d={45}^{\circ }$, $\phi =0^{\circ }$, ${\sigma }_{cl}$ = 20 km${}^{-1}$, and $r_{surf}$ = 0–1. Results of comparison of the calculated and interpolated data are shown in Figure 4. From Figure 4, it can be seen that the calculated and interpolated data coincide at the nodal points. The average difference between the results of calculations and interpolations in Figure 4a was $\Delta R$ = 2.9 km at $\lambda$ = 0.41 $\mathsf{\mu }$m, $\Delta R$ = 1.5 km at $\lambda$ = 0.55 $\mathsf{\mu }$m, and $\Delta R$ = 0.4 km at $\lambda$ = 0.86 $\mathsf{\mu }$m. The average differences between the calculated and interpolated results shown in Figure 4b was $\Delta R$ = 0.4 km at $\lambda$ = 0.41 $\mathsf{\mu }$m, $\Delta R$ = 0.1 km at $\lambda$ = 0.55 $\mathsf{\mu }$m, and $\Delta R$ = 0.05 km at $\lambda$ = 0.86 $\mathsf{\mu }$m. The reason for the differences between the calculated and interpolated values at intermediate grid nodes is the following. The values of $I_{sum,cloud}\left(R\right)$, $I_{sun,clear}$, ${\tilde{I}}_{surf,clear}$, $E_{0,clear}$, and ${\gamma }_{1,clear}$ were calculated by the Monte Carlo method with a statistical error of less than 0.01%. However, for $\Delta r_{surf}$ values close to 0.005, the derivative $d\Delta r_{surf}/dR$ is small; therefore, the computing error of the Monte Carlo method introduces a significant error in the definition of $R_{*}$. In particular, this error is the main reason for the observed discrepancy between the calculations and interpolation in Figure 4. The special feature in the behavior of $\Delta r_{surf}$ is clearly shown in Figure 5a. An analogous special feature in the behavior of $\Delta r_{surf}$ can also be seen in Figure 11 of the work [6] for distances from the cloud greater than 2–3 km.

From Figure 4b, it can be seen that with increasing $r_{surf}$, the $R_{*}$ values for some situations decrease, and for others increase. The reason is that, on the one hand, with increasing $r_{surf}$, the relative contribution of $\Delta r_{surf}$ = 0.005 to the $r_{surf}$ value decreases, and on the other hand, $I_{sum,cloud}$ increases with $r_{surf}$. On the whole, from Figure 4, it can be seen that the proposed interpolation formula solves well the problem of determining the $R_{*}$ value. The proposed approximation formula can be used for five MODIS channels (Nos. 1–4 and 8) under the following observation conditions: presence of cumulus cloudiness, cloud cover indices ${\delta }_{cl}$ from 0 to 0.5, ground surface reflectances $r_{surf}$ from 0 to 1, aerosol optical depths of the cloudless atmosphere ${\tau }_{0.55}$ from 0 to 0.89, cloud extinction coefficients in the range from 10 to 40 km${}^{-1}$, solar zenith angles and zenith angles of the receiving system from 0 to 60${}^{\circ }$, and azimuthal angle $\phi$ from 0 to 180${}^{\circ }$. Calculations were carried out for cloudless midlatitude summer LOWTRAN-7 models. As shown by the preliminary estimates, the difference between the results obtained can reach 10–20% for different models of the cloudless atmosphere, but the same AOD. In future, we plan to expand the possibilities of the algorithm for a wider set of models of cloudiness and cloudless atmosphere.

5. Approbation Method

For approbation of the interpolation formulas, two fragments of the MODIS image MOD021KM.A2017172.0325.006.2017172133827.hdf recorded on June 21, 2017 were analyzed. These two fragments were chosen because the average AOD of cloudless areas was significant for the first fragment, whereas for the second fragment, it was much smaller. Calculations were carried out for the fourth MODIS channel (centered at $\lambda$ = 0.55 $\mathsf{\mu }$m). The average values used in calculations of the fragments and the fragment boundaries are presented in

Table 1. Boundaries of the examined image fragments and average optical and geometrical conditions according to the MODIS data.

Parameter	Fragment 1	Fragment 2
Coordinates	53.4–56.4 N, 109–115 E	49.0–51.0 N, 121–123 E
Cloud cover index ${\delta }_{cl}$	0.15	0.087
Average altitude of the upper boundary of cloudiness $h_{max}$, km	4.1	2.6
Average optical thickness of cloudiness ${\tau }_{cl}$	30	15
Average AOD of cloudless areas	1.25	0.43
Average reflectance ${\overline{r}}_{surf}$ from the MODIS data	0.071	0.046
Average solar zenith angle ${\theta }_{sun}$, deg	34	27
Average satellite zenith angle ${\theta }_d$, deg	28	34
Average azimuthal angle $\phi$, deg	152	166

. The average cloud size was $\overline{L}$ = 1 km, and the lower cloud boundary was at the altitude $h_{min}$ = 1 km. In this case, $I_{sum,cloud}\left(R\right)$ entering into Formula (6) was calculated for P = 300 and N = 5000. For all results obtained, $\delta I\le$ 0.01.

The $\Delta r_{surf}\left(R\right)$ values calculated from Formula (7) are shown in Figure 5.

Based on the results shown in Figure 5 and the interpolation formula (12), the $R_{*}$ values were calculated. The results obtained are presented in

Table 2. Comparison of $R_{*}$ values calculated and interpolated for two test image fragments.

Image Fragment	Calculation, km	Interpolation, km	Difference, km
1	17	15	2
2	3.5	2	1.5

. The mask of pixels located at distances $R_{*}$ and less from the boundaries of the cloud projections onto the ground surface is shown in Figure 6. The comparison of the calculated and interpolated $R_{*}$ values differed by 2 km for fragment 1 and by 1.5 km for fragment 2. Thus, the interpolation allows one to estimate the radius of the cloud adjacency effect. Moreover, the computational time for the interpolation was, on average, 750 times less than with the application of the Monte Carlo method. This allows one to use the proposed method for problem solving in real time.

6. Conclusions

In this work, the method of estimating the cloud adjacency effect on the accuracy of the reconstructed reflectance of the ground surface areas non-shadowed by clouds based on calculation of the cloud field intensity averaged over an ensemble of realizations in the center of the entrance pupil of the receiving system has been proposed. The approbation of the method for two test fragments of the MODIS image MOD021KM.A2017172.0325.006. 2017172133827.hdf recorded in the fourth MODIS channel shows that the proposed interpolation formula allows the radius $R_{*}$ of the adjacency effect to be estimated. The difference between the values of the gap radius $R_{*}$ calculated by the Monte Carlo method and obtained by interpolation did not exceed 2 km for these areas.

In this stage of investigations, the interpolation formula for $R_{*}$ has been constructed for the midlatitude summer LOWTRAN-7 model and the OPAC model of cumulus cloudiness. In future, we plan to expand the possibilities of the algorithm to include other optical models of the atmosphere and cloudiness. The reliability of the results is provided by the performed test comparisons and the use of the well-known models of the atmosphere and cloudiness.

Verification of the results of this method is planned in the next stage of our research. It is suggested to use the following approach. For verification, a series of three images of the same cloudless fragments of the ground surface will be selected. Two images should not contain cloudiness at distances less than 50 km from the examined area, and in the third image, the cloudiness should be observed near the examined area. Then, analyzing the differences between the reconstructed reflectances, the distance from the cloudiness can be determined, for which these differences would be less than 0.005.

The advantages of the proposed approach over existing methods are that it allows one to estimate the adjacency effect averaged over an ensemble of cloud field realizations (cloud arrangement) on the reconstructed reflectance of cloudless ground surface areas with preliminary consideration of a wide range of optical and geometrical conditions of observations. Consideration of the deterministic gap allows more realistic situations to be examined compared to homogeneous cloud fields without deterministic gaps. In contrast to the above-mentioned papers [10,11], the proposed approach has no restrictions on the ground surface reflectance.