The Bidirectional Reflectance Factor (BRF) is defined as the ratio of the surface-reflected radiance to radiance reflected from an ideal Lambertian surface into the same beam geometry and illuminated by the same mono-directional beam (Martonchik et al., 2000; Schaepman-Strub et al., 2006). It describes the magnitude and angular distribution of surface reflected radiation in the absence of atmosphere and varies with the directions to the Sun, Ω 0 ∼ ( θ 0 , φ 0 ) , and to the sensor, Ω ∼ ( θ , φ ) . In this paper, the directions are expressed in terms of zenith, θ 0 and θ , and azimuthal, φ 0 and φ , angles. We will use symbols μ 0 and μ for cos ⁡ θ 0 and cos ⁡ θ , respectively.

For sufficiently dense vegetation such as equatorial forests, the BRF can be accurately approximated as (Knyazikhin et al., 2013) B R F λ ( Ω 0 , Ω ) = ρ ( Ω 0 , Ω ) i 0 1 − p × ω λ ( 1 − p ) 1 − p ω λ = D A S F ( Ω 0 , Ω ) × W λ   . (1)

The first factors on the right-hand side of Eq. 1 is the Directional Area Scattering Factor (DASF), which describes the canopy BRF if the foliage does not absorb radiation. The spectrally invariant DASF is a function of canopy geometrical properties, such as the tree crown shape and size, spatial distribution of trees on the ground, and within-crown foliage arrangement (Knyazikhin et al., 2013). The second factor, W λ , is the Canopy Scattering Coefficient (CSC), i.e., the fraction of intercepted radiation that has been reflected from, or diffusively transmitted through, the vegetation (Smolander and Stenberg 2005; Lewis and Disney 2007). The spectrally varying CSC is weakly sensitive to variation in the sun-sensor geometry. It conveys information about leaf optical properties (Knyazikhin et al., 2013; Latorre-Carmona et al., 2014; Adams et al., 2018).

Our forest BRF is parameterized in terms of spectrally invariant parameters (Knyazikhin et al., 2011; Stenberg et al., 2016). Here i 0 is the canopy interceptance defined as the portion of photons from the incident solar beam that collide with foliage elements for the first time. The symbol ρ designates the directional escape probability, i.e., the probability by which a photon scattered by a foliage element will exit the vegetation in the direction Ω through gaps. Spherical integration of π − 1 ρ | cos ⁡ θ | results in 1 − p , where p is the recollision probability, defined as the probability that a photon scattered by a foliage element in the canopy will interact within the canopy again. The spherical integration significantly weakens the sensitivity of p to sun-sensor geometry. Finally, ω λ is the wavelength dependent leaf albedo, i.e., the fraction of radiation incident on a leaf surface that is reflected or transmitted.

The directional escape probability controls the shape of the BRF. Indeed, photons scattered by sunlit leaves will escape the vegetation in the retro-illumination direction with unit probability since their paths are free of foliage elements. Photon paths in off-backscattering directions are more likely obstructed by leaves and the likelihood of photons escaping the canopy is consequently reduced. We follow methodology developed in (Yang et al., 2017) to simulate the hot spot effect. Kuusk’s model of the hot spot incorporated into the extinction coefficient of the radiative transfer equation is used to estimate the escape probability (Supplementary Appendix SA).

Our primary objective is to derive DASF from Terra MISR and DSCOVR EPIC observations. For vegetation canopies with a dark background, or sufficiently dense vegetation where the impact of canopy background is negligible, the DASF can be directly retrieved from the BRF spectrum in the weakly absorbing spectral interval, without involving canopy reflectance models, prior knowledge, or ancillary information regarding leaf scattering properties. We follow methodology developed in (Marshak and Knyazikhin 2017; Song et al., 2018) to approximate this variable using BRFs at NIR and green spectral bands: DASF is the ratio R / ( 1 − s ) where R and s are intercept and slope of the line passing two points ( B R F λ ω λ ,   B R F λ ) ,   λ =   green, NIR. Thus, D A S F ( Ω 0 , Ω ) = B R F green B R F NIR B R F green − β ( B R F NIR − B R F green )   .   (2)

Here β = ( 1 − ω NIR ) ω green / ( ω NIR − ω green ) where ω NIR and ω green represent leaf albedo of the brightest leaf at NIR and green spectral bands integrated over bandwidths. Its values are ω 555 = 0.461 , ω 865 = 0.978 ( β = 0.0196 ) for MISR and ω 551 = 0.490 ω 779 = 0.966 ( β = 0.035 ) for EPIC. These values were obtained from Lewis and Disney’s approximation (Lewis and Disney, 2007) of the PROSPECT model (Féret et al., 2008) with the following parameters: chlorophyll content of 16 μg cm⁻²; equivalent water thickness of 0.005 cm⁻¹, and dry matter content of 0.002 g cm⁻¹.

The probability of photons escaping the vegetation canopy depends on scattering order. The directional escape probability in Eq. 1 is an average over scattering orders (Supplementary Appendix SB). We approximate ρ ( Ω 0 , Ω ) by probabilities calculated for single scattered photons (Supplementary Appendix SC). We use the inclination index of foliage area to parameterize the leaf normal distribution (Ross 1981). This index characterizes the deviation of leaf orientation from the spherical distribution. It allows us to approximate the geometry factor, G, that appears in (Supplementary Appendix SA7 as G = 0.5 α , where the weight α varies between 0 and 2. The leaf normals, Ω L , are simulated by spherical distribution corrected for the deviation, i.e., g ( Ω L ) = α . The corresponding scattering anisotropy (Supplementary Appendix SA2) becomes: j ( Ω ) = 1 3 π ( sin ⁡ ϑ − γ ⁡ cos ⁡ ϑ ) + τ L 3 cos ⁡ ϑ 0.5 π   , (3) where ϑ = π − acosΩΩ 0 is the scattering angle (the angle between incident and scattered radiation) and τ L represents the leaf transmittance, which was set to 0.5 in our calculations. Under these assumptions DASF in the upward directions rearranges to the form (Supplementary Appendix SC) D A S F ≈ ρ 1 i 0 1 − p ( 1 ) = 1 1 − p ( 1 )   j ( Ω ) μ 0 μ   1 − t ψ L ψ   .   (4)

Here p ( 1 ) is the single scattering approximation of the recollision probability (Supplementary Appendix SC); t = exp ( − 0.5 ) = 0.61 ; ψ = μ 0 − 1 + μ − 1 ( 1 − ϰ H S ) ; the factor ϰ H S is defined by Supplementary Appendix SA4, and L is an effective extinction coefficient. Thus, our model depends on two parameters. They are the hot spot parameter h that appears in ϰ H S and the effective extinction coefficient L. The former determines the shape of DASF, while the latter controls its magnitude.

Our study is focused on equatorial evergreen broadleaf forests that include Amazonian central rainforests (0°–10°S and 70°–60°W), Congo rainforests in Central Africa (5°S–5°N and 20°–30°E) and Southeast Asian rainforests (19.80°–26.57°N and 92.5°–105°E). Figure 1 shows locations of our study area. The seasonal transition between wet and dry seasons is a distinct feature of tropical rainforests, which leads to intra-annual patterns of leaf flushing and abscission.

About 95% of our Amazonian central rainforest is covered with terra firme rainforests (Nepstad et al., 1994). The average annual rainfall during the 2000–2019 period is about 2,600 mm. The seasonal cycle consists of a short dry season, June to October, and a long wet season thereafter.

The equatorial rainforests of Central Africa are the second largest and least disturbed of the biodiversly-rich and highly productive rainforests on Earth (Cook et al., 2020). Our study area includes central and part of western and northeast Congolian lowland forests. The Congo basin exhibits bimodal precipitation pattern and has two wet and two dry seasons per year (Yang et al., 2015). The wet seasons occur in March-April-May and September-October-November, while dry season months are December-January-February and June-July-August. The average annual rainfall over the past 2 decades is about 1761 mm.

Our third region consists of two sub-regions depicted as Region 3.1 (23.50–26.57°N and 92.5–98.62°E) and 3.2 (19.80–21.54°N and 97.93–105°E). The first one is a subtropical moist broadleaf forest ecoregion in Mizoram–Manipur–Kachin rain forests. It occupies the lower hillsides of the mountainous border region joining India, Bangladesh, and Burma (Myanmar). The average annual rainfall over the past 20 years is about 1,545 mm. The dry season is from October to April, and wet season is May to September. Region 3.2 represents a subtropical moist broadleaf forest ecoregion in Northern Indochina. The wet seasons occur in May to September while dry season months are October to April.

Various variables from several independent satellite sensors over our study area were used in this research. These include land cover maps and leaf area index (LAI) from the MODerate resolution Imaging Spectroradiometer (MODIS), precipitation from Tropical Rainfall Measuring Mission (TRMM), surface bidirectional reflectance factor (BRF) from Multi-angle Imaging SpectroRadiometer (MISR) on the Terra platform and BRF from Earth Polychromatic Imaging Camera (EPIC) on Deep Space Climate Observatory (DSCOVR).

MODIS land cover dataset. Collection 6 Terra and Aqua MODIS land cover product from 2001 to 2019 at yearly temporal frequency and 0.05° spatial resolution (Friedl and Sulla-Menashe 2015) was used to identify our study area. This product provides several classification schemes. The map of LAI classification scheme was adopted in this research. Figure 2 illustrates LAI classification scheme used by DSCOVR EPIC operational algorithm for the generation of Vegetation Earth System Data Record (WWW-VESDR 2021).

MODIS LAI datasets. Collection 6 Terra and Aqua MODIS LAI products (Myneni et al., 2015a; Myneni et al., 2015b) for the period February 2000 to December 2019 were used in this study. The LAI dataset provides 8-days composite LAI at 500-m spatial resolution. The C6 MODIS LAI product was evaluated against ground-based measurements of LAI and through inter-comparisons with other satellite LAI products (Yan et al., 2016; Yan et al., 2016).

TRMM precipitation dataset. Monthly precipitation data from the TRMM (3B43 version 7) at 0.25° spatial resolution for the period January 2000 to December 2019 (WWW-TRMM 2011) was used in this study. This dataset provides the best-estimate precipitation rate and root-mean-square precipitation-error estimates by combining four independent precipitation fields (Huffman et al., 2007).

DSCOVR EPIC MAIAC dataset. Level 2 DSCOVR EPIC Multi-Angle Implementation of Atmospheric Correction (MAIAC, version 1) surface BRF and aerosol optical depth (AOD) at 551 nm from 2016 to 2019 were also used. The EPIC instrument has provided imageries in near backscattering directions with the phase angle between 4° and 12° at ten ultra-violet to near infrared (NIR) narrow spectral bands until June 27, 2019, when the spacecraft was placed in an extended safe hold due to degradation of the inertial navigation unit (gyros). DSCOVR returned to full operations on March 2, 2020 after the navigation problem had been resolved. After March 2020 the range of phase has substantially increased towards backscattering reaching 2° (Lyapustin et al., 2021; Marshak et al., 2021).

The MAIAC BRF are available at four spectral bands; they are 433 (band width 3.0) nm, 551 (3.0) nm, 680 (2.0) nm and 780 (2.0) nm. Data are projected on a 10-km SIN grid and available at 65–110 min temporal frequency (WWW-MAIAC 2018). EPIC sees Amazonian rainforests between 11 UTC and 18 UTC, Congo forests between 5 UTC and 14 UTC and Southeast Asian rainforests between 1 and 7 UTC.

MISR datasets. The MISR sensor views each 1.1 km ground pixel symmetrically about the nadir in the forward and aftward directions along the spacecraft’s flight track. Image data are acquired with nominal view zenith angles relative to the surface reference ellipsoid of 0.0⁰ (camera An), 26.1⁰ (Af and Aa), 46.5⁰ (Bf and Ba), 60.0⁰ (Cf and Ca) and 70.5⁰ (Df and Da) in four spectral bands centered at 446 (band width 41.9 nm), 558 (28.6) nm, 672 (21.9) nm, and 866 (39.7) nm. MISR obtains global coverage between ±82° latitudes in 9 days (Diner et al., 1998; Diner et al., 1999). Level 2 version 3 MISR land surface (WWW-MISR_SURFACE 1999) and aerosol (WWW-MISR_AEROSOL 1999) products for the period of January 2016 to December 2019 over our study area were used. The surface reflectance parameter BRF in 9 view angles and four MISR spectral bands is at 1.1 km spatial resolution. The aerosol optical depth is available at 4.4 km spatial resolution. Both parameters are projected on Space Oblique Mercator (SOM) projection, in which the reference meridian nominally follows the spacecraft ground track.

Directions from ground pixel to MISR cameras form view lines on the polar plane, which are characterized by slope, k = tan ⁡ υ K , and intercept, b (Figure 2A). The slope is aligned with ground track and is roughly constant with υ K ≈ 15.5 ° . The intercept is associated with location of pixel within the MISR 360 km swath. We parameterize MISR BRF in terms of the solar zenith angle, phase angle and intercept. The phase angle, γ , is calculated as γ = acosΩΩ 0 = acos ( cos ⁡ θ ⁡ cos ⁡ θ 0 + sin ⁡ θ ⁡ sin ⁡ θ 0 ⁡ cos ( φ − φ 0 ) )   , (5) where Ω 0 ∼ ( θ 0 , φ 0 ) and Ω ∼ ( θ , φ ) are directions from ground pixel to the Sun and sensor, respectively. We assign the sign “plus” to the phase angle if the MISR view direction approaches to the direction to the Sun from North, i.e., sin ⁡ θ ⁡ cos ( υ K − φ ) > sin ⁡ θ 0 ⁡ cos ( υ K − φ 0 ) , and “minus” otherwise (Figure 2B).

The MODIS LAI and TRMM precipitation data over forested pixels were selected using flags indicating highest retrieval quality. The 8-days 500 m LAI products over our study area (Figure 1) were spatially aggregated to 0.01° and 0.1° resolutions which were then used in our analyses.

The MISR and DSCOVR EPIC surface BRF over our study area were first refined by removing pixels with aerosol optical depth over 0.3. MISR and EPIC datasets were further re-projected to 0.01° and 0.1° Climate Modeling Grids (CMG), respectively. For each pixel, MISR and EPIC DASFs were calculated using Eq. 2, which then were used to generate monthly DASFs. If there were several observations of a pixel within a given month, a median DASF value was assigned to such pixel.

Area-averaged DASF as a function of mean SZA and phase angle, γ , is defined as D A S F ( S Z A , γ ) = ∑ x y ∈ A D A S F x y ( S Z A x y , V Z A x y ) cos ⁡ S Z A x y ∑ x y ∈ A ⁡ cos ⁡ S Z A x y   , (6) where the summation is over pixels ( x , y ) in the selected area A at which the phase angle takes a given value γ .

Equation 4 is used to simulate DASF. It depends on the hot spot parameter, h, and effective extinction coefficient, L. The former is a function of SZA and determines the angular shape of DASF, while the latter controls its magnitude and depends on LAI. The following two-step fitting technique was implemented to derive equations for h and L using monthly MISR DASF.

Step 1: Matching angular shapes of observed and modeled DASF. For a given month, we used S Z A x y and monthly average MODIS LAI as a first approximation to the effective extinction coefficient (i.e., L x y ≈ L A I x y ) to simulate D A S F x y at each pixel ( x , y ) in MISR view angles as a function of h. Next, we used Eq. 6 to calculate area-averaged simulated-DASF as a function of hot spot parameter, h. Finally, we selected h that minimized ( R 2 − 1 ) 2 + ( s − 1 ) 2 , where R 2 and s are the coefficient of determination and slope of the relationship between area averaged values of observed and simulated DASFs. The selected hot spot parameter provides the best agreement between angular shapes of modeled and observed DASF.

Step 2: Matching magnitudes of observed and modeled DASF. For a given month, we used SZA and h (SZA) to simulate D A S F x y at MISR view angles as a function of L. Eq. 6 was used to calculate area averaged simulated-DASF as a function of effective extinction coefficient, L, i.e., D A S F ¯ m o d ( γ , L ) . We selected L that minimizes Normalized Root Mean Square Error (NRMSE) between simulated, D A S F ¯ m o d , and observed, D A S F ¯ M I S R , area-averaged DASFs, i.e., N R M S E ( L ) = 1 N ∑ γ [ D A S F ¯ m o d ( γ , L ) − D A S F ¯ M I S R ( γ ) ] 2 1 N ∑ γ D A S F ¯ M I S R ( γ , L )   → m i n .   (7) This value of L matches magnitudes of observed and modelled DASFs.

Monthly MISR DASF data for the 2017 to 2019 period over our study area (Figure 1) were used to execute our two-step fitting procedure. The SZA exhibits small variation within our regions during a month and therefore can be accurately represented by its monthly mean. A time series of the solutions to the Step-1 procedure therefore gives a set of the hot spot parameters corresponding to different SZA. Seasonal variations of LAI in equatorial forests allowed us to accumulate solutions to the Step-2 procedure corresponding to different values of MODIS LAI. We used those sets to derive dependences of the hot spot parameter and effective extinction coefficient on SZA and LAI, respectively.

Figure 3 shows an example of our two-step fitting technique for Congo forests (region 2) in September-2018. As illustrated in Figure 4, Eq. 4 approximates observed DASF to within NRMSE = 8% and R² = 0.85. The largest difference between observed and simulated DASFs occurred at phase angles above 90⁰. Such points are separated by an ellipse in Figure 4. For PA > 90⁰, MISR BRFs were mainly acquired by off-nadir F and D cameras, which have higher uncertainties compared to near nadir observations.

The sets of solutions to the Steps 1 and 2 procedures allowed us to regress the hot spot parameter, h, and effective extinction coefficient, L, versus SZA and MODIS LAI, respectively, as (Figure 5) h ( S Z A ) = 5.96 − 5.90 ⁡ cos ⁡ S Z A ≈ 11.92 ⁡ sin 2 S Z A 2 , (8) L = 2.93 ⋅ L A I M O D I S − 8.53 (9) There was no correlation between L and SZA, as expected.

Thus, our model for DASF of equatorial forests is generated by Eq. 4 with the hot spot parameter h and effective extinction coefficient L given by Eqs 8, 9. It has two input parameters; they are Sun position in the sky, Ω 0 ∼ ( θ 0 , φ 0 ) , and MODIS LAI.

Observed versus modeled DASF. We used monthly MISR DASF for the period between 2017 and 2019 to derive equations for the hot spot parameter and effective extinction coefficient. The proximity between observed and modeled DASFs were characterized by NRMSE = 8%, R² = 0.85 (Figure 4). We analyzed modelled and observed monthly DASF for Year 2016 to see if the performance metrics is similar to that of the training data set. Figure 6 illustrates monthly MISR DASF and its simulated counterpart for Amazonian forests in April 2016. The largest differences between them are at high phase angles. Figure 7 shows MISR DASF plotted versus modeled DASF accumulated over our study area during Year 2016. The comparison suggests a good performance of Eq. 4 to simulate MISR DASF over equatorial forests.

Diurnal variations of observed and modeled DASFs. The next step in the assessment of our approach is to see if the model can reproduce diurnal variation of monthly EPIC DASF. Figure 8 shows examples of diurnal variations in observed and modeled DASFs for 3 regions in our study area. As one can see the largest deviation between model and observation occurs when SZA exceeds 60⁰. The uncertainty of the MAIAC BRF product is low for the EPIC observations near the local noon. It however may significantly increase at high zenith angles resulting in an underestimation of surface BRF (Lyapustin et al., 2021). And this is what we see in Figure 8.

Scatter plot of diurnal courses of modeled and EPIC DASFs accumulated over Amazonian, Congo forests and region 3.2 in southeast Asia during the 2017 to 2019 period is shown in Figure 9. Note that data from December to November over Congo forests are not present in this plot. For these regions, our model approximates diurnal courses of the observed DASFs to within NRMSE = 7% with R² = 0.82.

On average, modeled DASF over Congo during December through February overestimates observed DASF by about 20%. About 70% of data on the scatter plane are located within a 15% circle centered at mean values of EPIC and modeled DASFs and therefore differ from respective mean values by less than 15%.

For the region 3.1 in Southeast Asian rainforest, modeled DASF overestimates observations by about 5%. The data are also concentrated on the scatter plane: about 75% of data on the model-vs.-observation scatter plane are concentrated within a 15% circle centered at mean values of observed and modeled DASFs. The R² is consequently low (Y = 0.8X+0.08, R² = 0.32).

In summary, Eq. 4 can accurately reproduce DASF in terms of proximity to both angular variations observed by MISR and diurnal courses measured by DSCOVR EPIC sensor. It therefore provides a strong basis for synergy of DSCOVR EPIC and Terra MISR sensors to monitor changes in equatorial forests. Our next step is to see if the model can detect changes.

The forest structural organization determines the magnitude and angular variation of DASF (Schull et al., 2011; Knyazikhin et al., 2013). Its angular signatures therefore provide unique and rich sources of diagnostic information about forests. Here we analyze DASF over our study area to see if it can detect seasonal changes of the equatorial forests.

The seasonal transition between wet and dry seasons is a distinct feature of equatorial rainforests, which leads to intra-annual patterns of leaf flushing and abscission (Samanta et al., 2012; Bi et al., 2015). Since our study is focused on structurally intact and undisturbed regions of the equatorial forests (i.e., no changes in forest geometry), variation in leaf area is a key factor causing variation in DASF.

We start with analyses of variation in the DASF acquired over Amazonian central rainforest. In situ studies and satellite data indicated higher leaf area during the dry season relative to the wet season (Huete et al., 2006; Hutyra et al., 2007; Myneni et al., 2007; Hilker et al., 2014; Jones et al., 2014). The growth-limiting impact of water deficit on rainforest during the dry season is alleviated through deep roots and hydraulic redistribution (Oliveira et al., 2005; Pierret et al., 2016), resulting in a sunlight mediated seasonality in leaf area (Bi et al., 2015). Figure 10 illustrates these findings, that is, green leaf area increases during the dry season (June to October), has high values during the early part of the wet season (November to October) and decreases thereafter (March to May).

Let us compare observed DASFs from the late dry season (October) and middle part of the wet season (March). Eq. 4 predicts that an increase in the effective extinction coefficient, with SZA unchanged, increases the magnitude of DASF at all phase angles, i.e., results in an upward shift in the angular signature of the DASF, as illustrated in Figure 3. The SZAs in the select region of Amazonian forests in March (SZA = 25.5, std = 1.2) and October (SZA = 20.5, std = 1.1) are very close. At low SZA such a small difference minimally impacts the shape of angular signatures. As one can see in Figure 11, both MISR and EPIC show a distinct decrease in DASF in all phase angles between October and March with no discernible change in the overall shape of the angular signatures. Such a simple change in magnitude can only result from a change in LAI since other structural variables, such as tree crown shape and size do not vary seasonally in this forest.

Let us now consider DASF in the early (June) and late (October) dry seasons. LAI has changed from about 5.5 to 6.4. MISR and EPIC measurements are made at significantly higher SZA in June (SZA = 37.6, std = 2.3) compared to October (SZA = 20.6, std = 1.1). The magnitude and shape of angular signatures are impacted when both canopy properties and SZA vary as middle panel in Figure 11 illustrates. This makes the comparison of the signatures difficult. We can transform the June’s signature to the sun-sensor geometry in October using Eq. 4. As right panel of Figure 11 demonstrates the transformed DASF is a downward shift of the October’s DASF, indicating a lower LAI in June. DASFs of the remining regions in our study area exhibit similar behavior (not shown here).

Our next example demonstrates time series of DSCOVR EPIC DASF acquired over the Congolese forests. The Congo basin exhibits bimodal precipitation pattern and has two wet and two dry seasons per year (Yang et al., 2015). The wet seasons occur in March -April-May and September-October-November, while dry season months are December to February and June to August. Unlike Amazonian forests, monthly average LAI follow the patterns of precipitation (Figure 12). It exhibits notable bimodal seasonal variations.

The above analyses have demonstrated that an increase in LAI, with SZA unchanged, results in an upward shift in the angular signature of the DASF (Figure 3). The EPIC DASF at fixed solar zenith and phase angles therefore should covary with LAI. The Earth-observing geometry of the EPIC sensor is characterized by phase angle between 2⁰ and 12⁰. A question then arises whether or not such small variation in phase angle can be ignored. Therefore, we examine two algorithms to generate EPIC time series. The first one selects EPIC observations at SZA = 25⁰, 30⁰ and 45⁰ irrespective of values of the phase angle. If there are no reflectance data under these illumination conditions during a month, we transform DASF to a desired SZA. In the second case, we select observation at fixed sun-sensor geometries. Figure 13 shows LAI and DASF at fixed SZA = 30⁰ and varying phase angle. At low SZA, the EPIC time series correlates well with the bimodal seasonal variation of LAI, as expected. This also suggests that Eq. 4 is an effective tool to fill missing data at a given fixed SZA.

An increase in SZA however can eliminate the bimodal feature of DASF. This is illustrated in Figure 14 showing annual courses of EPIC DASF generated by the two algorithms introduced above. As one can see in left panel of this figure, the EPIC time series at SZA = 45⁰ becomes flat between May and October. Two factors are responsible for this effect. First, the decrease of phase angle to its local minimum in July enhances DASF and therefore tends to suppress decrease in DASF due to the dry season decrease in LAI. At low SZA LAI has a stronger impact on DASF than phase angle. The impact of phase angle however increases with SZA and can become a dominant factor causing variation in DASF. In our example, this occurs at a SZA of 45⁰ and higher. As right panel of Figure 14 illustrates, DASF at fixed SZA and phase angle retains its bimodal property. Thus, both SZA and phase angle should be taken into account when analyzing DSCOVR EPIC data. Eq. 4 therefore becomes of particular importance for analyses of EPIC observations over vegetated land. A strong effect of phase angle on EPIC reflectance was recently documented in (Marshak et al., 2021). Our analyses reinforce this effect.