Vegetation Angular Signatures of Equatorial Forests From DSCOVR EPIC and Terra MISR Observations

In vegetation canopies cross-shading between finite dimensional leaves leads to a peak in reflectance in the retro-illumination direction. This effect is called the hot spot in optical remote sensing. The hotspot region in reflectance of vegetated surfaces represents the most information-rich directions in the angular distribution of canopy reflected radiation. This paper presents a new approach for generating hot spot signatures of equatorial forests from synergistic analyses of multiangle observations from the Multiangle Imaging SpectroRadiometer (MISR) on Terra platform and near backscattering reflectance data from the Earth Polychromatic Imaging Camera (EPIC) onboard NOAA’s Deep Space Climate Observatory (DSCOVR). A canopy radiation model parameterized in terms of canopy spectral invariants underlies the theoretical basis for joining Terra MISR and DSCOVR EPIC data. The proposed model can accurately reproduce both MISR angular signatures acquired at 10:30 local solar time and diurnal courses of EPIC reflectance (NRMSE < 9%, R2 > 0.8). Analyses of time series of the hot spot signature suggest its ability to unambiguously detect seasonal changes of equatorial forests.


INTRODUCTION
The global forest ecosystem absorbs about 25% of the total anthropogenic CO 2 emission from atmosphere via carbon accumulation to forest biomass (Reichstein et al., 2013). Forests store 75% of terrestrial carbon, and account for 40% of the carbon exchange with atmosphere each year (Schlesinger and Bernhardt, 2012). Within the forest ecosystem, tropical forests contain about 40-50% of the terrestrial carbon stock (Lewis et al., 2009) and are potentially responsible for about 70% of terrestrial carbon sink (Pan et al., 2011). Monitoring and quantifying changes in tropical forests therefore play a critical role in understanding the global carbon cycle and future climate change.
Monitoring of dense vegetation such as equatorial rainforests represents the most complicated case in optical remote sensing because reflection of solar radiation saturates and becomes weakly sensitive to vegetation changes. At the same time, the satellite data are strongly influenced by changing sun-sensor geometry. This makes it difficult to discriminate between vegetation changes and sun-sensor geometry effects. For instance, studies on Amazon forest seasonality based on analyses of data from single-viewing sensors disagree on whether there is more greenness in the dry season than in the wet season: the observed increase in vegetation indices were explained by an increase in leaf area, an artifact of sun-sensor-geometry and changes in leaf age through the leaf flush (Huete et al., 2006;Brando et al., 2010;Samanta et al., 2012;Morton et al., 2014;Saleska et al., 2016). The impact of droughts on Amazon forests has also been debated Samanta et al., 2010;Samanta et al., 2011;Xu et al., 2011). Conflicting conclusions among these studies arose from different interpretations of surface reflectance data acquired under saturation conditions (Bi et al., 2015). Developing methodologies that allow us to unambiguously interpret reflectance of dense forests is worthy of special attention.
Broadly used approaches for interpretation of satellite data from single-viewing sensors consider the viewing and solar zenith angle dependence of reflected radiation to be a problematic source of noise or error, requiring a correction or normalization to a "standard" sun-sensor geometry (Lyapustin et al., 2018). Transformation of such data to a fixed standard sunsensor geometry therefore invokes statistical assumptions that may not apply to specific scenes. The lack of information about angular variation of forest reflected radiation introduces model uncertainties that in turn may have significant impact on interpretation of satellite data (Gorkavyi et al., 2021).
Unlike single-angle methodologies, multiangle approaches exploit angular variation of surface reflected radiation as unique and rich sources of diagnostic information and enable the rigorous use of the radiative transfer theory. In vegetation canopies cross-shading between finite dimensional leaves leads to a peak in reflectance in the retro-illumination direction. This effect is called the hot spot in optical remote sensing (Gerstl and Simmer, 1986;Ross and Marshak, 1988;Kuusk, 1991;Myneni, 1991). The hotspot region in reflectance of vegetated surfaces represents the most information-rich directions in the angular distribution of canopy reflected radiation. The hot spot phenomenon correlates with canopy architectural parameters such as foliage size and shape, crown geometry and withincrown foliage arrangement, foliage grouping, leaf area index and its sunlit fraction (Ross and Marshak, 1991;Qin et al., 1996;Goel et al., 1997;Qin et al., 2002;Yang et al., 2017;Pisek et al., 2021). Angular signatures that include the hot spot region are critical for monitoring phenological changes in equatorial forests (Bi et al., 2015). Availability of hot spot signatures of equatorial forests would make monitoring their changes more reliable.
The Multiangle Imaging SpectroRadiometer (MISR) on Terra platform provides simultaneous multiangle observations of surface reflectance since December 1999. Its observing strategy allows for a good angular variation of surface reflectance in equatorial zone. However spatially and temporally varying phase angle 1 could be far from zero, making frequent observations of canopy reflectance in the hot spot region impossible. The NASA's Earth Polychromatic Imaging Camera (EPIC) onboard NOAA's Deep Space Climate Observatory (DSCOVR) was launched on February 11, 2015 to the Sun-Earth Lagrangian L1 point where it began to collect radiance data of the entire sunlit Earth every 65-110 min in June 2015. It provides imageries in near backscattering directions .
The DSCOVR EPIC observations therefore provide unique information required to extend angular sampling of the MISR sensor to the hot spot region. The objectives of this paper are to 1) develop a new methodology that synergistically incorporates features of Terra MISR and DSCOVR EPIC observation geometries and results in hot spot signatures of equatorial forests; 2) generate angular signatures of equatorial rainforests for the period of concurrent Terra MISR and DSCOVR EPIC data and asses their quality; 3) demonstrate their value for monitoring seasonal changes of the equatorial forests.

Reflectance of Dense Vegetation
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) (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 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 1 The phase angle is the angle between the directions to the Sun and sensor Frontiers in Remote Sensing | www.frontiersin.org October 2021 | Volume 2 | Article 766805 2 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 sunsensor 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 BRF λ ω λ , BRF λ , λ green, NIR. Thus, 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 −2 ; equivalent water thickness of 0.005 cm −1 , and dry matter content of 0.002 g cm −1 .

Approximation of DASF
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: 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) Here p (1) is the single scattering approximation of the recollision probability (Supplementary Appendix SC); t exp(−0.5) 0.61; ψ μ −1 0 + μ −1 (1 − 8 HS ); the factor 8 HS 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 8 HS and the effective extinction coefficient L. The former determines the shape of DASF, while the latter controls its magnitude.

Study Area
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.  and phase angle (PA), the latter is the angle between the directions to the Sun and sensor. We assign the sign "plus" to the phase angle if the MISR view direction approaches the direction to the sun from North (i.e., above the dotted red line perpendicular to the MISR view line), and "minus" otherwise.
Frontiers in Remote Sensing | www.frontiersin.org October 2021 | Volume 2 | Article 766805 4 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.

Data Used
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 intercomparisons with other satellite LAI products .
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 bestestimate precipitation rate and root-mean-square precipitationerror 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 0 (camera An), 26.1 0 (Af and Aa), 46.5 0 (Bf and Ba), 60.0 0 (Cf and Ca) and 70.5 0 (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, c, is calculated as c acosΩΩ 0 acos cos θ cos θ 0 + sin θ sin θ 0 cos φ − φ 0 , 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).

Data Processing
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 Frontiers in Remote Sensing | www.frontiersin.org October 2021 | Volume 2 | Article 766805 5 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, c, is defined as DASF SZA, c xy∈A DASF xy SZA xy , VZA xy cos SZA xy xy∈A cos SZA xy , where the summation is over pixels (x, y) in the selected area A at which the phase angle takes a given value c.

Hot Spot Parameter and Effective Extinction Coefficient
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 SZA xy and monthly average MODIS LAI as a first approximation to the effective extinction coefficient (i.e., L xy ≈ LAI xy ) to simulate DASF xy 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 DASF xy 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., DASF mod (c, L). We selected L that minimizes Normalized Root Mean Square Error (NRMSE) between simulated, DASF mod , and observed, DASF MISR , area-averaged DASFs, i.e., 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 2 0.85. The largest difference between observed and simulated DASFs occurred at phase angles above 90 0 . Such points are separated by an ellipse in Figure 4. For PA > 90 0 , 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(SZA) 5.96 − 5.90 cos SZA ≈ 11.92 sin 2 SZA 2 , L 2.93 · LAI MODIS − 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.

Assessment of DASF
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 2 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 0 . The uncertainty of the MAIAC BRF product is Frontiers in Remote Sensing | www.frontiersin.org October 2021 | Volume 2 | Article 766805 6 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 2 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 2 is consequently low (Y 0.8X+0.08, R 2 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.

Monitoring Equatorial Forests
The forest structural organization determines the magnitude and angular variation of DASF 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 intraannual 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; FIGURE 3 | (A) MISR DASF of region 2 (Congo forests) in September-2018 (hollow circles). Its step-1 and step-2 approximations are shown as crosses and dots, respectively. The dashed line is a polynomial fit to the Step-2 approximation. (B) MISR DASF versus step-1 (crosses) and step-2 (circles) approximations. NRMSEs between MISR DASF and its Step 1 and Step 2 approximations are 12 and 4%, respectively. Mean SZA (std) 21.2 0 (1.5 0 ), h 0.8, LAI 5.6, L 7.15.     (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 0 and 12 0 . 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 0 , 30 0 and 45 0 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 0 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 0 becomes flat between May and October. Two factors FIGURE 9 | Correlation between 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. Data from December to February over Congo forests are excluded. NRMSE 7%. The relationship between these data is characterized by a regression line with a slope of 1 and negligible intercept; R 2 0.82. 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 0 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 . Our analyses reinforce this effect.

SUMMARY AND CONCLUSIONS
We used Directional Area Scattering Function (DASF) to characterize angular signatures of equatorial forests. It describes the canopy BRF if the foliage does not absorb radiation and is a purely structural variable. For vegetation canopies with a dark background, or sufficiently dense vegetation where the impact of canopy background is negligible, the DASF can be accurately approximated from the BRF in the weakly absorbing spectral intervals without involving canopy reflectance models, prior knowledge, or ancillary information regarding leaf scattering properties (Knyazikhin et al., 2013). Equation 2 is used to obtain approximations of DASF from the Terra MISR and DSCOVR EPIC data. The DASF becomes independent on spectral band composition of a sensor acquiring surface reflectance data, which is an important prerequisite for achieving consistency and complementarity between DSCOVR EPIC and Terra MISR observations. We adapted a model for the canopy hot spot implemented in the operational algorithm for generation of Earth System Data Record (VESDR) from DSCOVR EPIC observations (Yang et al., 2017;WWW-VESDR 2021). In this approach, the sunlit leaves are treated as a stochastic reflecting boundary, which depends on distribution of leaves in the canopy space and the Sun position in the sky. Photons reflected by the boundary can either enter the vegetation canopy or exit it. The shaded leaves represent the interior points. Their interactions with photons are described by a stochastic radiative transfer equation. The directional escape probability that appears in   Eq. 1 is a weighted sum of photons reflected by the boundary and canopy interior points. Kuusk's model of the hot spot incorporated into the extinction coefficient (Supplementary Appendix SA) is used to evaluate the escape probability as a function of scattering order, which is then used to calculate the average escape probability (Supplementary Appendix SB).
Contributions of multiple scattered photons are accounted by the recollision probability.
Here we simplified this model. First, a one-dimensional radiative transfer equation is used to simulate canopy radiative regime (Supplementary Appendix SA6). Second, the average escape probability is approximated by a probability calculated for first order scattered photons (Supplementary Appendix SC). Under these assumptions, DASF is approximated by a simple equation that depends on two parameters. They are the hot spot parameter that appears in the canopy hot spot coefficient and the effective extinction coefficient. The former determines the shape of DASF, while the latter controls its magnitude. These two parameters should be specified to generate angular signatures of equatorial forests.
In spite of substantial theoretical advancement in modeling the radiative transfer in vegetation canopies, quantitative data on the hot spot are still few and far between. Here we specified the hot spot parameter by fitting shapes of observed and modeled DASF using MODIS LAI as an initial approximation to the effective extinction coefficient. The hot spot parameter was found to be almost proportional to 1 − cos SZA (R 2 0.96) with a coefficient of proportionality around 6 (left panel in Figure 5). The trigonometric term can be interpreted as a correction of the canopy hot spot coefficient (Supplementary Appendix SA4) for errors due to its approximation by a constant value (Supplementary Appendix SA) whereas the coefficient of proportionality as a mean linear dimension of foliage elements (Knyazikhin and Marshak 1991;Nilson 1991) specific to equatorial forests. This equation for the hotspot parameter was used in all our calculations.
The effective extinction coefficient determines the magnitude of the DASF. Theoretically this variable can be obtained by replacing the 3D extinction coefficient with an effective value that provides a best agreement between horizontally averaged canopy reflectances and solutions of 1D radiative transfer equations. Basically, it depends on LAI, leaf normal distributions and clumping indices. In our approach, the effective extinction coefficient was obtained by matching magnitudes of MISR and shape-adjusted modeled DASFs. This coefficient was found to be linearly related to the MODIS LAI (R 2 0.65, right panel in Figure 5). We used this relationship in all our calculations.
Note the MODIS LAI was used as a first approximation to the effective extinction coefficient, which then was iterated to its optimal value. Alternatively, one can use relationships between LAI and various vegetation induces (e.g., NDVI) to make rough estimates of LAI first and then iterate them to the extinction coefficient. This procedure may result in relationships between the extinction coefficient and vegetation indices, which can make the model dependent on the hot spot parameter and vegetation indices.
Our model for angular signatures of equatorial forests can accurately reproduce both MISR angular signatures acquired at 10:30 local solar time and diurnal course of EPIC reflectance (NRMSE<9%, R 2 > 0.8) and therefore assures consistency and complementarity between DSCOVR EPIC and Terra MISR observations. This provides a powerful tool to argue for changes in vegetation structure as it was demonstrated in our analyses of seasonal variations of angular signatures acquired over equatorial forests.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: DSCOVR EPIC and Terra MISR data were obtained from the NASA Langley Research Center Atmospheric Science Data Center (ASDC) (https://epic.gsfc.nasa.gov). The MODIS products were acquired from the Land Processes (LP) Distributed Active Archive Center (DAAC) (https://lpdaac.usgs. gov). TRMM 3B43 data are publicly available from the NASA