Separating overstory and understory leaf area indices for global needleleaf and deciduous broadleaf forests by fusion of MODIS and MISR data

Forest overstory and understory layers differ in carbon and water cycle regimes and phenology, as well as ecosystem functions. Separate retrievals of leaf area index (LAI) for these two layers would help to improve modeling forest biogeochemical cycles, evaluating forest ecosystem functions and also remote sensing of forest canopies by inversion of canopy reflectance models. In this paper, overstory and understory LAI values were estimated separately for global needleleaf and deciduous broadleaf forests by fusing MISR and MODIS observations. Monthly forest understory LAI was retrieved from the forest understory reflectivity estimated using MISR data. After correcting for the background contribution using monthly mean forest understory reflectivities, the forest overstory LAI was estimated from MODIS observations. The results demonstrate that the largest extent of forest understory vegetation is present in the boreal forest zones at northern latitudes. Significant seasonal variations occur for understory vegetation in these zones with LAI values up to 2–3 from June to August. The mean proportion of understory LAI to total LAI is greater than 30 %. Higher understory LAI values are found in needleleaf forests (with a mean value of 1.06 for evergreen needleleaf forests and 1.04 for deciduous needleleaf forests) than in deciduous broadleaf forests (0.96) due to the more clumped foliage and easier penetration of light to the forest floor in needleleaf forests. Spatially and seasonally variable forest understory reflectivity helps to account for the effects of the forest background on LAI retrieval while compared with constant forest background. The retrieved forest overstory and understory LAI values were compared with an existing dataset for larch forests in eastern Siberia (40–75 N, 45–180 E). The retrieved overstory and understory LAI is close to that of the existing dataset, with an absolute error of 0.02 (0.06), relative error of 1.3 % (14.3 %) and RMSE of 0.93 (0.29) for overstory (understory). The comparisons between our results and field measurements in eight forest sites show that the R2 values are 0.52 and 0.62, and the RMSEs are 1.36 and 0.62 for overstory and understory LAI, respectively.


Introduction
Forests not only provide habitat and food for animals and fiber and fuel for human beings but also control the global climate and biogeochemical cycles.Most forests, excluding tropical rainforests, have two distinct layers: an overstory layer that mainly consists of arbors, and an understory layer that includes shrubs, grasses and moss.These layers are often treated differently in ecosystem modeling due to their different photosynthetic capacity, carbon residence times, phenology and climatic and environmental responses (Vogel and Gower, 1998;Rentch et al., 2003;Marques and Oliveira, 2004;Kim et al., 2016).The understory vegetation plays an essential role in supporting biodiversity, the nutrient cycling and the capability of soil and water conservation Published by Copernicus Publications on behalf of the European Geosciences Union.
Y. Liu et al.: Separating overstory and understory leaf area indices in forests (Suchar and Crookston, 2010;Qiao et al., 2014).The leaf area index (LAI) is critical for describing the water, carbon and energy exchange of vegetation with the atmosphere (Braswell et al., 1997;Gitelson and Kaufman, 1998).A global wall-to-wall LAI dataset with separation of forest LAI for overstory and understory layers would help to improve the modeling of forest carbon and water cycles and the evaluation of forest ecosystem functions (Law and Waring, 1994).
Satellite remote sensing provides powerful tools for the estimation of global forest LAI.Several global LAI products have been produced from satellite observations, such as MODIS LAI and fractional photosynthetically active radiation product MOD15 from Terra-Aqua/MODIS data (Myneni et al., 2002), CYCLOPES (Baret et al., 2007) and GEOV-1 (Baret et al., 2013) from SPOT/VEGETATION data, MERIS LAI from ENVISAT/MERIS data (Bacour et al., 2006) and GLOBMAP global long-term LAI from a combination of Terra/MODIS and NOAA/AVHRR data (Liu et al., 2012a).However, the LAI in these products is the total LAI, i.e., the sum of the forest overstory and understory LAI.Considerable efforts have been made to separate forest overstory and understory LAI at a regional scale.For example, the overstory and understory LAI values were estimated by using a combination of high spectral and high spatial resolution images using a neural network method and field measurements in the Longmenhe forest nature reserve in China (Huang et al., 2011).LAI was also estimated separately for larch forests in northern Eurasia using its relationship with the normalized difference water index (NDWI) based on the three-dimensional radiative transfer model (Kobayashi et al., 2010), which assumed that the understory LAI is stable over the entire year; this means that it is only applicable to deciduous forests with evergreen understory.As forests can vary considerably even within a particular land cover type due to their complex structures and multiplicity of species, a global dataset is highly desirable to describe the heterogeneous spatial distribution of forest overstory and understory.
Multi-angle remote sensing can capture signals of different forest layers because the observed proportions for different forest layers vary with the viewing angle, making it possible to separate forest overstory and understory LAI on a global scale.Forest background reflectivity, which is the reflectance of materials below the forest canopy, including understory vegetation, leaf litter, moss, lichen, rock, soil, snow, etc., was estimated from the Multi-angle Imaging Spectroradiometer (MISR) (Canisius and Chen, 2007;Pisek and Chen, 2009) and two-angle observations with Compact Airborne Spectrographic Imager (Pisek et al., 2010a) observations based on the four-scale model (Chen and Leblanc, 1997).If the forest background reflectivity was known, LAI could be estimated separately for the forest overstory and understory.The forest background reflectivity over North America has been derived from MISR observations at a spatial resolution of 1 • (Pisek and Chen, 2009), and that result has been used to correct the effects of forest background for the determination of the forest overstory LAI over North America (Pisek et al., 2010b).
Background reflectivities in the red and near-infrared (NIR) bands have been retrieved over global forest areas at a resolution of 1.1 km using MISR observations from 2000 to 2010 (Jiao et al., 2014).A large amount of invalid retrievals exists in these background reflectivity maps for each day because of frequently missing data and invalid retrievals in MISR land surface products, which is due to cloud and aerosol contamination, topographically complex terrain, and other suboptimal illumination conditions.To generate spatially complete 1 km resolution maps, the monthly mean forest background reflectivity was derived by averaging all available valid retrievals for each day during the 11-year period.This dataset makes it possible to separate the forest overstory and understory LAI globally.
The forest overstory is the main component for carbon fixation and changes seasonally and interannually.If the forest background condition is assumed to change only from month to month but remain stable over the period of different years, the high temporal resolution of overstory LAI can be inferred from satellite measurements with the help of the MISR forest background reflectivity.It takes 9 days for the MISR to acquire global coverage, which makes it challenging to capture the integrated seasonal patterns of global forest overstory.As MODIS provides global radiative measurements within 1 day, the combination of MODIS data and MISR forest background reflectivity could help estimation of forest overstory LAI.In this paper, the forest overstory and understory LAI values were separated for global needleleaf and deciduous broadleaf forest based on the global 1 km MISR forest background reflectivity and MODIS observations from 2008 to 2010.The monthly mean forest understory LAI was retrieved from the MISR forest background reflectivities in the red and NIR bands.The overstory LAI was determined from the MODIS land surface reflectance with the effects of forest background accounted for by using the MISR forest background reflectivity.Because notable uncertainties may be introduced in the understory LAI retrieval for evergreen broadleaf forests (EBFs) in tropical zones due to the excessive complexity of the forest structure, which makes it difficult to separate it into two layers, the evergreen broadleaf forest type was excluded in this study.

MODIS and MISR background reflectivity data
The forest overstory LAI was retrieved from the MODIS land surface reflectance product MOD09A1 from 2008 to 2010.MOD09A1 provides global 8-day composite 500 m resolution land surface reflectance in bands 1-7 without the effects of the atmospheric gases and aerosols since February 2000.The land cover type was defined by the MODIS land cover type product MCD12Q1, which supplies a yearly land cover classification map with a 500 m resolution of the globe derived through a supervised decision tree classification method with five different classification schemes.In this study, the product from the International Geosphere-Biosphere Programme (IGBP) global vegetation classification scheme is selected.The MOD09A1 and MCD12Q1 products are all provided on the sinusoidal grid.
The forest understory LAI was derived from the MISR forest background reflectivity maps by Jiao et al. (2014).This MISR seasonal forest background reflectivity dataset provides monthly reflectivities in the red and NIR bands over global forest areas at 0.01 • (approximately 1 km) resolution on a geographic grid, which represents the reflectance of all materials below the forest canopy, such as understory vegetation, rocks, soils, leaf litter, lichens, mosses, snow or their mixture.The reflectivity was derived from the MISR surface bidirectional reflectance factor (BRF) in the nadir and 45.6 • forward directions (An and Bf cameras) for each day from 2000 to 2010 based on the four-scale model.Then, monthly forest background reflectivity was produced by averaging 11year valid results for each month to replace a large number of invalid retrievals (Jiao et al., 2014).
All MODIS data were re-projected to a geographic grid, which is the same as the forest background reflectivity.MOD09A1 and MCD12Q1 images from 2008 to 2010 were transformed to a geographic reference system at a 0.005 • × 0.005 • (approximately 500 m) spatial resolution using nearest neighbor interpolation and then composed to form global maps.The MOD09A1 land surface reflectance was screened for cloud contamination based on a refined cloud mask for MODIS land surface reflectance products, which identifies cloudy pixels based on the inflexion point between clear-sky and cloudy observations of time series of reflectances assemblage for the same location (Liu and Liu, 2013).The snow/ice pixels were also labeled with MOD09A1 state flags.

Field measurements of forest overstory and understory
Ground LAI measurements of forest overstory and understory were collected from relevant published papers for the evaluation of the derived LAI products.These include 28 field LAI measurements of forest overstory at eight sites and 12 measurements of forest understory at Prince Albert National Park, Saskatchewan, Canada.Detailed information about these measurements is presented in Table 1.The overstory LAI measurements provided effective LAI in five conifer-dominated boreal forest sites in Finland, while true LAI was provided for other two deciduous broadleaf forest and one deciduous needleleaf forest sites along a wide latitu-dinal gradient in the Northern Hemisphere (37.75-66.45• N).
Here, those effective LAI values were converted to true LAI using the clumping index of the corresponding pixel in the global 500 m resolution map derived from MODIS bidirectional reflectance distribution function (BRDF) product (He et al., 2012).Furthermore, field measurements of spectral hemispherical-directional reflectance factors of forest understory at seven sites were transferred to understory normalized difference vegetation index (NDVI), and also used to evaluate the retrieved understory LAI.These sites include different forest types and understory species: a sparse black spruce forest in Alaska (PFRR), a dense black spruce forests in Canada (Sudbury), two southern boreal forest stands in Finland with dominant tree species of Scots pine (Hyytiälä xeric) and birches (Hyytiälä herb rich), hemiboreal needleleaf (pine, Järvselja radiation transfer model intercomparison (RAMI) pine) and deciduous (birch, Järvselja RAMI birch) stands in Estonia, and a temperate mixed forest in Switzerland (Laegern).Table 2 shows the detailed information of these sites.The understory vegetation is mainly composed of herbaceous species at the Hyytiälä herb rich, Järvselja RAMI birch and Laegern sites, while it is composed of shrubs, herbs and mosses at the other four sites.

Estimation of LAI
In this study, the forest overstory refers to the tree canopy, while the forest understory refers to the vegetation below the forest overstory tree canopy, mainly including shrubs, grasses and moss.Forest overstory and understory LAI was estimated using an algorithm for the separation of forest LAI between overstory and understory layers based on the LAI algorithm of the GLOBCARBON project of European Space Agency (ESA) (Deng et al., 2006).The derived forest overstory and understory LAI datasets are named GLOBMAP forest overstory LAI (hereafter referred to as LAI o ) and GLOBMAP forest understory LAI (LAI u ), respectively.The GLOBCARBON LAI algorithm produces the LAI using the land-cover-dependent relationships between the LAI and the vegetation index (VI) with consideration of the BRDF effects explicitly based on the four-scale model and Chebyshev polynomials (Deng et al., 2006).Forests are combined to four biomes (coniferous, tropical, deciduous and mixed), and non-forest vegetation is combined into two biomes, i.e., shrub and grass/crop/other vegetated surfaces.Different algorithm coefficients are utilized for each biome type based on separate four-scale simulations.The relationship based on the reduced simple ratio (RSR) (Brown et al., 2000) is used for forests, while the relationship based on the simple ratio (SR) (Jordan, 1969) is used for shrub, grasses and other nonforest biomes.First, the effective LAI (L E ) is derived based on the function of VI (SR or RSR): where f L E _VI is the biome-specific function defining the relationships between L E and BRDF-modified VI at a specific view and sun angle combination (θ v , θ s , ϕ); here, SR is used for non-forest biomes and RSR for forest biomes.VI obs refers to satellite-observed VI (SR or RSR).f biome is the function defining the different algorithms for forest (f forest ), shrub (f shrub ) and grass (f grass ) biomes, which are presented in Sect.2.2.1 and 2.2.2.Function f BRDF , which quantifies the BRDF effects of VI obs , depends on the angular reflectance behavior of the spectral bands involved.Next, the true LAI is calculated from effective LAI retrievals (L E ) based on clumping index ( ), which accounts for the vegetation clumping effect on the plant and canopy scales: (2) In this paper, the uses the global clumping index map derived from the MODIS BRDF product (He et al., 2012).The forest understory LAI was estimated from the MISR forest background reflectivity based on the LAI algorithms for grass and shrub (see Sect. 2.2.1).The forest overstory LAI was derived from the MODIS land surface reflectance data with the effects from the background corrected based on the MISR monthly forest background reflectivity.Then, the forest total LAI was calculated by summing the forest overstory LAI and understory LAI. Figure 1 shows the general flowchart for the forest overstory and understory LAI separation algorithm.The LAI derived from the MISR forest background reflectivity (Jiao et al., 2014) should approximate the forest understory LAI.Here, by assuming that the forest understory has a similar composition to a mixture of shrubland, grassland and moss, its LAI was retrieved from the SR of the forest background based on the LAI algorithms for shrub and grass, which is applicable for grass and other non-forest biomes except shrub.First, the vegetation index SR for the forest background (SR B ) was calculated using MISR background reflectivity in the red (ρ Red_B ) and NIR bands (ρ NIR_B ): Then, the forest understory effective LAI was retrieved based on the function of SR B using the LAI algorithm for shrub and grass, respectively (Eq. 1).The VI uses SR, and the function f L E _VI in Eq. ( 1) uses the relationship between L E and SR for shrub and grass.The effective LAI is directly derived from the SR (here using SR B ) with no corrections except for BRDF (Eqs. 4 and 5).The corresponding solar zenith angle (θ s ), view zenith angle (θ v ) and relative azimuth angle (ϕ) of the MISR cameras, which were used to generate the MISR forest background reflectivity, were also used for correcting the BRDF effects in the retrieval of the forest understory LAI (Eq.1).Because the observations from two MISR cameras were used in the generation of the MISR forest background reflectivity, including the nadir and 45.6 • forward cameras, the forest understory effective LAI was calculated using the geometries in these two cameras separately, and the average value of the two results was used.
After that, the true LAI was calculated from effective LAI retrievals using the global mean value of clumping index for shrubs (0.73) and grasses (0.75) that were derived from MODIS BRDF products (He et al., 2012) (Eq.2).Then, these two retrievals for the two biomes were averaged and used as the true LAI for the forest understory (LAI u ).Finally, the global monthly LAI u maps were produced by averaging the valid LAI u from 2000 to 2010 for each month to generate spatially complete maps.The monthly LAI u values are missing only 0.07 % of global needleleaf and deciduous broadleaf forests areas due to missing MISR background reflectivity during 2000-2010.In addition, some of derived LAI u values exceed the valid range (0-6).These invalid retrievals occurred mainly in summer, which are probably due to the large uncertainties in background reflectivity for dense canopy.These missing values and invalid retrievals were flagged in LAI u maps.

Estimation of the forest overstory LAI
The effects of the forest background on the MODIS observations were accounted for with the monthly forest background reflectivity maps.Next, the forest overstory LAI was retrieved using modified MODIS observations.In the GLOB-CARBON LAI algorithm, the effects of the background are considered by using the background SR value of 2.4 in all of the simulations of the four-scale model for all forest types (Deng et al., 2006).To account for the effect of the difference between the standard background SR value (2.4) used in the model simulation and the actual background SR (SR B ), which varies among the sites and seasons, the MODIS observed SR was adjusted by removing the effects of the background (Pisek et al., 2010b): where SR Overstory_modified is the adjusted SR for the forest overstory; SR obs refers to observed SR, which is calculated from the MODIS land surface reflectance (MOD09A1) in the red and NIR bands; and SR MAX is the maximum SR value of the algorithm for a forest type at the view zenith angle (θ v ).
The monthly forest background SR (SR B ) was calculated from the monthly MISR forest background reflectivity in the red and NIR bands.As there are still missing values in monthly SR B maps, and uncertainties exist in MISR forest background reflectivity dataset, especially for dense canopy (Jiao et al., 2014), the monthly SR B was scaled to a 10 km resolution to fill the missing values and reduce these uncertainties and then used to adjust the MODIS observed SR (SR obs ) based on Eq. ( 6).Since SR B represents the signals from all materials below the forest canopy, the effects of the forest background were accounted for.After that, the IGBP land forest classes in the MODIS land cover product (MCD12Q1) were grouped into four forest biomes (conifer, tropical, deciduous and mixed forest).The overstory effective LAI was retrieved from this adjusted MODIS SR (SR Overstory_modified ) and corresponding MODIS geometry data using the LAI algorithm for each forest type (Eq.1).Since RSR is less sensitive to variable background of forest stands (Brown et al., 2000), the function f L E _VI in Eq. ( 1) uses the relationship between L E and RSR (Eq.7) for global evergreen needleleaf forests (ENFs), deciduous needleleaf forests (DNFs) and deciduous broadleaf forests (DBFs), and the effective LAI was retrieved from RSR with consideration of the BRDF effects in the shortwave infrared band (SWIR, MODIS band 5) (Eq.8).where ρ SWIR is the SWIR (band 5) reflectance for MODIS; the function f SWIR_BRDF , quantifying the BRDF effects of MODIS band 5, depends on the angular reflectance behavior of this band; and ρ SWIRmax and ρ SWIRmin are the maximum and minimum values of the SWIR reflectance.Finally, the true LAI for the forest overstory (LAI o ) was calculated from the effective overstory LAI based on Eq. ( 2) using the global clumping index map derived from the MODIS BRDF product (He et al., 2012).

Calculation of the forest total LAI
The forest total LAI (GLOBMAP forest total LAI, hereafter referred to as LAI T ) was calculated by summing the forest overstory LAI and the understory LAI from the same month: (9)

Estimation of the GLOBCARBON LAI
The LAI was also retrieved from the MODIS land surface reflectance products MOD09A1 for the global forest areas based on the GLOBCARBON LAI algorithm (hereafter referred to as GLOBCARBON LAI) to evaluate the LAI o , LAI u and LAI T .In the GLOBCARBON LAI algorithm, the constant background SR value of 2.4 is used for all forest types without consideration of the spatial and seasonal variations of the forest background (Deng et al., 2006).This constant background SR value is based on field measurements from boreal forests in Canada, which includes effects due to green mosses and the understory (Chen et al., 2002).It does not differentiate forest LAI between overstory and understory.The GLOBCARBON LAI is approximately equal to the forest total LAI with some of the effects due to the forest understory considered (Deng et al., 2006).3 Results

Sensitivity analysis of forest understory vegetation composition
Since it is challenging to acquire the composition of understory vegetation at a global scale, the understory LAI was calculated by averaging the LAI retrievals from MISR background reflectivity based on the GLOBCARBON LAI algorithm for non-forest biomes, i.e., shrub and grass/crop/other vegetated surfaces (see details in Sect.2.2.1).In fact, the structure of understory vegetation is usually complex and spatially disparate in contrast to half shrubs and half other non-forest vegetation (Peltoniemi et al., 2005).Thus, the smaller the difference of the understory LAI retrieved based on the algorithm between shrub (LAI u,shrub ) and grass (LAI u,grass ), the less sensitive the LAI u is to the composition of understory vegetation, which may result in less uncertainties in derived LAI u .In this section, the difference between LAI u,shrub (Eq.4) and LAI u,grass (Eq.5) was calculated to evaluate the sensitivity of derived LAI u to forest understory vegetation composition for global ENFs, DNFs and DBFs.
Figure 2 shows the distribution of differences between LAI u,grass and LAI u,shrub for global ENFs, DNFs, DBFs and all three types of forests combined.The two LAI results are close to each other with the difference concentrated around 0 www.biogeosciences.net/14/1093/2017/Biogeosciences, 14, 1093-1110, 2017 and 82 % (94 %) of values within ±0.2 (±0.3).LAI u,grass is slightly greater than LAI u,shrub for most forest pixels.The annual mean difference is 0.11 for global forests except EBFs, and DBFs show slightly larger difference (0.12) than ENFs (0.09) and DNFs (0.11).Seasonal cycles were present in the difference between LAI u,grass and LAI u,shrub .Figure 3a and b show temporal profiles of monthly mean difference between these two algorithms for various forest types in the Northern and Southern hemispheres, respectively.Generally, the difference is larger in summer.In the Northern Hemisphere, the difference is less than 0.10 due to the sparse understory vegetation from October to March, while it reaches the maximum values during May and June with mean difference ranging from 0.20 to 0.28.The three forest types show a noticeable discrepancy in the seasonal curve of mean differences.A large difference is concentrated in June in DNFs, with the largest monthly mean difference up to 0.28.In the Southern Hemisphere, the mean difference is generally smaller than in the Northern Hemisphere.For ENFs and DNFs, the monthly mean difference is generally below 0.10.For DBFs, the mean difference shows significant seasonal variations with the opposite seasonal curve compared to the Northern Hemisphere.
A large mean difference is present from August to March, with most of the monthly mean difference ranging from 0.12 to 0.19.The maximum monthly mean difference appears in December (0.19).This suggests that the uncertainties in derived LAI u are generally slightly larger in summer and in the Northern Hemisphere, especially for DNF in the Northern Hemisphere in June.

Global distribution of forest overstory and understory LAI
The monthly LAI u was estimated from the MISR monthly background reflectivity over the global needleleaf and DBF area at a spatial resolution of 1 km, and the 8-day LAI o was retrieved from the MODIS land surface reflectance data MOD09A1 from 2008 to 2010 at a spatial resolution of 500 m.
Figure 4 shows the global distribution of LAI u for each month, and Fig. 5 presents global LAI o maps for each month from 2010.Significant spatial and seasonal variations exist for forest understory and overstory LAI.Forest understory is mostly found in the boreal forest zone in the northern latitudes (50-70 • N), especially in summer.In July, 84 % of valid retrievals are concentrated in this zone, and this percentage is still up to 54 % in January, when much of the understory vegetation disappears or is covered by snow.The LAI u is sparse for the southern latitudes (23.5-63.0• S), where the forest area represents only approximately 4 % of global forests.Only 2 % of the valid LAI u retrievals are found in this region in July, and this percentage is less than 7 % even in January.The LAI u values mainly range from 0 to 3, with most values occurring in the range 0-1.These values could be up to 2-3 in July and August in boreal forests.Notable seasonal variations are found for LAI u , especially for boreal forests at high latitudes in the Northern Hemisphere.In this region, the LAI u is mainly less than 0.5 from November to April.Its value increases starting in May and reaches approximately 1.0 in May with a maximum value of up to 2-3 from June to August.Then, the values decrease in September and are below 1.0 in October.
The spatial and seasonal patterns of the LAI o are close to that of the forest total LAI, suggesting that overstory is the dominant component of forest LAI (Fig. 5).Compared with the southern latitudes, the seasonal variations in LAI o are much more pronounced in the region of 30-70 • N, where deciduous forests are widely distributed.The LAI o is approximately 0 from November to March due to snow covers, and its maximum value is greater than 4.0 from June to August.In the southern latitudes (23.5-63.0• S), the valid retrievals of the LAI o are sparse with inconspicuous seasonal variations due to the small forest area.

Distribution of overstory and understory LAI for different forest biomes
The LAI o and LAI u values were compared for three major forest types, including ENFs, DNFs and DBFs.The histograms and mean values of monthly LAI u and 8-day LAI o (2008)(2009)(2010) were analyzed over the global forest area for each forest type.Figure 6a and b show histograms of the LAI o and LAI u , respectively.For the forest overstory, the LAI values range from 0 to 10, with the majority concentrated in the range of 0-8.Broadleaf forests have more LAI o values greater than 4 than needleleaf forests.The percentage with LAI o above 3 is up to 29.3 % for DBFs, while this percentage is 24.1 and 18.0 % for ENFs and DNFs, respectively.The mean value of the LAI o is also larger for DBFs (2.17) than for ENFs (2.12) or DNFs (1.72).This result is mainly attributed to its short www.biogeosciences.net/14/1093/2017/Biogeosciences, 14, 1093-1110, 2017 growing season and the low solar radiation at high latitudes, where DNFs are mainly distributed.
For the forest understory, the LAI values are smaller than those for the overstory; most LAI u values are in the range from 0.0 to 3.0.In contrast to LAI o , the needleleaf forests more LAI u values greater than 1.5 than the DBFs; 19 % of LAI u values are above 1.5 for DBFs, while this percentage is 24 % for ENFs and 23 % for DNFs.The mean value of the global LAI u is also larger for needleleaf forests (1.06 for ENFs and 1.04 for DNFs) than for DBFs (0.96).This result is probably attributed to the more clumped foliage (Chen et al., 1997) and generally sparse tree spacing of the needleleaf forests, leading to more radiation penetrating through the canopy than for the less-clumped broadleaf forests, which is favorable for the growth of understory vegetation (Gower et al., 1999).In addition, the signals from the understory vegetation may be relatively easier to observe by satellite sensors for highly clumped needleleaf forests.

Seasonal variations of the forest overstory and understory LAI
The temporal profiles for LAI o and LAI u for the DBFs, ENFs and DNFs in the Northern Hemisphere are presented in this section to further evaluate these products.For each forest type, LAI u values of all pixels covered with a specific forest type were averaged for each month.Similarly, LAI o was averaged for each 8-day interval during the period of 2008-2010.These mean LAI values were calculated for the Northern Hemisphere, and the seasonal curves are presented in this section, where the majority of forests in the world are found.
Figure 7 shows time series of the overstory (thick lines) and understory (thin lines) LAI.Significant seasonal variations of the LAI o and LAI u are presented for all three forest types, and the shapes of the curves are similar between the LAI o and LAI u for each type.The differences in the magnitude of the seasonal variations and the lengths of the growing seasons are presented for three types.For the LAI o , DBFs show a larger magnitude of seasonal variation and a longer growing season than needleleaf forests, with LAI o values approximately 0.6 from November to March and up to more than 4.5 in June and July.For ENFs, the LAI o values in the summer (approximately 3.0) are smaller than those of broadleaf forests.Although the LAI o values for ENFs are larger than those of deciduous forests in the winter, the values are still below 2.0, which is probably due to the low sensitivity of satellite measurements in the red and NIR bands to the canopy LAI as a result of the decrease in the foliage chlorophyll in the winter.The magnitude of the seasonal variations is also the smallest for DNFs compared to the other two forest types, where the LAI o is only approximately 3.0 in July.Because DNFs are usually found at high latitudes or in mountains with low temperatures and precipitations, the lengths of the growing seasons for DNFs are generally shorter than those of DBFs and ENFs.The seasonal curves for LAI u are similar to those for LAI o with smaller magnitudes.In contrast to the forest overstory, the needleleaf forest understory shows higher LAI values (up to 2.2 for DNFs and 1.9 for ENFs in July) than that for broadleaf forests (1.7 for DBFs) in the summer, and thus presents larger seasonal variation than that of broadleaf forests, especially for DNFs.The mean proportion of understory LAI to total LAI is generally greater than 30 %, indicating that the understory cannot be ignored.This percentage is 36 % for DBFs, which is slightly larger than that for DNFs (34 %) and ENFs (31 %).Differences in the lengths of the growing seasons for LAI u are also found for various forest types, which is similar to the forest overstory.

Comparison with the overstory and understory LAI of larch forests in northern Eurasia
The derived LAI o and LAI u were compared with an existing overstory and understory LAI dataset, which was estimated for the larch forest with understory vegetation mainly composed of evergreen shrubs, such as cowberry and other deciduous species, covering eastern Siberia (40-75 • N, 45-180 • E) (Kobayashi et al., 2010).This understory LAI dataset (Kobayashi LAI u ) was derived through the relationship with the NDWI on the day of leaf appearance.The overstory LAI dataset (Kobayashi LAI o ) was estimated from the relationships between the overstory LAI and the seasonal in-  creases in the NDWI after leaf appearance.The algorithm was applied to SPOT/VGT observations and produced understory and overstory LAI for the larch forests defined in Global Land Cover database for the 2000 (GLC2000) deciduous needleleaf class over eastern Siberia with the resolutions of 1/112  8 shows the density scatter plots of the two LAI datasets, and the statistical analysis results are shown in Table 3.For understory LAI, both datasets are concentrated in LAI ranges of 0.0-1.0.The comparison shows good agreement with a high density of points dispersed along the 1 : 1 line and a root mean square error (RMSE) of 0.29.The mean GLOBMAP LAI u (0.48) is also close to that of Kobayashi LAI u (0.42), leading to an absolute error of 0.06 and a relative error of 14.3 %.For overstory LAI, the majority of the points are distributed in the range from 0.0 to 4.0, which represents the notable seasonal variations of the foliage in these DNFs.The area with the highest density (dark red) is also along the 1 : 1 line, especially in the LAI range from 0.0 to 3.0, which indicates that the majority of the GLOBMAP LAI o agrees with the Kobayashi LAI o .The mean value of the GLOBMAP LAI o is 1.95, and this value is 1.61 for the Kobayashi LAI o , resulting in an absolute error of 0.34 and a relative error of 21.1 %.Relatively high densities are also found in the Kobayashi LAI o of 4.0, which can probably be attributed to the relatively concentrated distributions of high LAI values in maximum values in the Kobayashi LAI algorithms.Some points with relatively low densities (blue) stray from the 1 : 1 line, and the mean value of the GLOBMAP LAI o is slightly higher than that of the Kobayashi dataset, indicating that the Kobayashi LAI o is slightly smaller than the GLOBMAP LAI o .The maximum value of the Kobayashi LAI o is 4.0, based on the statistics from 2008 to 2009, while the maximum value of the GLOBMAP LAI o is 10.0, with the majority in the range from 0.0 to 6.0.To avoid the influence from this difference, a statistical analysis was performed using only the pixels with overstory LAI ranging 0.0 to 4.0 for the two datasets.The mean values of the overstory LAI for Y. Liu et al.: Separating overstory and understory leaf area indices the GLOBMAP LAI o (1.53) and the Kobayashi LAI o (1.51) are very similar, with an absolute error of 0.02 and a relative error of 1.3 %, and decreased RMSE from 1.42 to 0.93.

Validation with field measurements
It is challenging to directly validate the satellite LAI products with spatial resolutions from several hundreds of meters to kilometers against ground measurements due to the uncertainties from scaling, heterogeneity, geolocation and the limited spatial and temporal sampling of ground data (Privette et al., 2000;Weiss et al., 2007).When it comes to the forest overstory and understory LAI, this work is even more challenging due to the lack of separate field measurements for these two layers, especially for understory LAI.In this section, the GLOBMAP LAI u and LAI o were compared with field LAI measurements extracted from related references at eight forest sites that cover major forest types except for EBFs.Additionally, the understory LAI (estimated field LAI u ) values were also estimated from the field understory NDVI measurements (field NDVI u ) at the seven sites (Table 2), and then compared with the GLOBMAP LAI u .The algorithms for shrub and for grass/crop/other vegetated surfaces were applied, and the mean value of the two retrievals was used as understory LAI.To reduce the uncertainties from heterogeneity and geolocation, the LAI retrievals on 3 × 3 km pixels around the site were averaged and then compared with the ground data.For those field measurements that provide mean LAI for a month or a period of time, the mean values of LAI retrievals during these periods were calculated for comparison.
For forest understory, the LAI measurements (N = 7) at Prince Albert National Park, Canada (53.70 • N, 106.20 • W), in April to October from 2000 to 2003 were extracted from figures in Barr et al. (2007).The mean values of field understory LAI were calculated for each month from April to October and then compared with the monthly mean GLOBMAP LAI u at corresponding geolocations.Figure 9a shows the results of comparison.The plots are close to the 1 : 1 line, with a slope of 0.96 and an offset of 0.17.The R 2 is 0.62 and RMSE is 0.62.For forest overstory, the LAI measurements (N = 28) at eight sites from 2000 to 2013 were compared with the GLOBMAP LAI o derived from MOD09A1 land surface reflectance obtained in the nearest time periods.Figure 9b shows the comparison results for forest overstory LAI.GLOBMAP LAI o seems to slightly overestimate the ground data, with a slope of 1.22 and an offset of −0.07.The R 2 is 0.52 and RMSE is 1.36.The LAI o generally agrees well with the ground data for overstory LAI less than 2.0, with plots close to the 1 : 1 line.On the contrary, the differences between the LAI o and ground data increase when overstory LAI is greater than 2.0.Similarly, the differences between the LAI u and the field measurement are also smaller for understory LAI values less than 1.5 than those that are approximately 2.0 (Fig. 9a).These results could be probably attributed to the larger uncertainties in reflectivities of the forest background for dense canopies (Jiao et al., 2014), which will affect the understory and overstory LAI retrieved for forest with high LAI values.In addition, the differences in the temporal period for field measurements and satellite observation also affect the comparisons, especially for understory LAI which uses the multi-year average values.
Figure 10 presents the time series of GLOBMAP LAI u , field NDVI u and estimated field LAI u at seven forest sites.Generally, GLOBMAP LAI u shows reasonable seasonal curves and is consistent well with the shapes of field NDVI u and estimated field LAI u .Notable seasonal variations are represented with GLOBMAP LAI u , which could characterize the seasonal curves of dominated deciduous understory vegetation at these sites except at Laegern.At the northernmost boreal forest site, PFRR (65.12 • N, 147.5 • W), the length of growing season for understory is short, with GLOBMAP LAI u greater than 1.0 only in June, July and August.The GLOBMAP LAI u overestimates the estimated field LAI u at this site.At the most southern site, Sudbury (47.16 • N, 81.76 • W), although the overstorydominated species (black spruce) is the same as at PFRR, the GLOBMAP LAI u is greater than 1.0 until October.Our result in July is very close to the estimated field LAI u in day of year (DOY) 177.
Hyytiälä herb rich (61.84 • N, 24.32 • E) and Hyytiälä xeric (61.81 • N, 24.33 • E) are two southern boreal forest stations with dominant tree species of birches and Scots pine, respectively.The former is more fertile than the latter.Understory vegetation is dominated by herbaceous species and graminoids at the herb-rich site, while lichens and heather dominate the xeric heath forest.At the Hyytiälä herb rich site, the GLOBMAP LAI u captures the seasonal curves represented by field NDVI u , and the values are consistent well with estimated field LAI u , with the LAI u reaching approximately 2.0 in the summer.In contrast, at the Hyytiälä xeric site, the LAI u shows smaller variations with the values below 1.5 in the summer.
For the two hemiboreal forests sites, Järvselja RAMI birch (58.28 • N, 27.33 • E) and RAMI pine (58.31 • N, 27.30 • E), the geolocations are more southern, and the understory vegetation is more abundant and hosts more species than in boreal Hyytiälä.The GLOBMAP LAI u increases earlier and the values are still above 0.5 in October in the more southern Järvselja sites.As the complexity of the understory enhances the challenge of understory LAI retrieval, the seasonal series are not as smooth as that in Hyytiälä.Even so, the GLOBMAP LAI u could capture the understory development in the spring and senescence in the fall, and its values are generally close to the estimated field LAI u .
The Laegern site in Switzerland represents a temperate mixed forest.The overstory is dominated by very big European beech and Norway spruce, with effective LAI up to 5.5, mean tree height of 30.6 m and maximum crown radius greater than 10 m (see details in Pisek et al., 2016).The un-  derstory vegetation is sparse and consists mainly of Allium ursinum.Although the GLOBMAP LAI u in September (0.9) is very close to the estimated field LAI u in DOY250, its time series does not appear stable during the whole year.This suggests that large uncertainties may exist in the retrieved understory and overstory LAI for such dense and closed canopy which make it difficult to observe the signals from understory on the top of canopy, which has also been pointed out by Pisek et al. (2016).

Seasonal effects of the background reflectivity on the LAI retrieval
The GLOBMAP LAI T was generated by summing LAI o and LAI u in the same month, and the mean differences between the GLOBCARBON LAI and the GLOBMAP LAI T and LAI o (GLOBCARBON LAI minus GLOBMAP LAI T /LAI o ) were calculated over the global forest area from 2008 to 2010 to evaluate the effects of monthly pixel-specific forest background reflectivity on forest LAI retrieval.
Figure 11 shows maps of the difference between the GLOBCARBON LAI and the GLOBMAP LAI T as well as LAI o over the global forest area.The difference between the GLOBCARBON LAI and the GLOBMAP LAI T is negative for most forest pixels.This is probably because some forest understory effects have been considered in the GLOBCAR-BON LAI by using the constant background SR value of 2.4 (Deng et al., 2006), while the GLOBMAP LAI T includes the LAI for all forest canopy and understory vegetation.This difference is small for needleleaf forests, with a yearly mean difference of −0.48 for ENFs and −0.59 for DNFs, while the difference is larger for DBFs (−0.68).This suggests that the standard background SR in the GLOBCARBON LAI algorithm is more suitable for boreal forests, which is probably because the standard background SR value (2.4) is based on field measurements in boreal forests in Canada (Chen et al., 2002).In contrast, the difference between the GLOBCAR-BON LAI and the GLOBMAP LAI o is positive for most forest pixels, which is attributed to the partial signals of understory still included in the GLOBCARBON LAI.Because the difference mainly represents signals from the forest understory, it is also larger for needleleaf forests than broadleaf forests due to their more clumped foliage, especially in the boreal forest zones (Fig. 11b), with yearly mean differences of 0.62 for ENFs, 0.43 for DNFs and 0.35 for DBFs. Figure 12 shows monthly mean differences and standard deviation (SD) series in the Northern Hemisphere.Notable seasonal variations are also found for the three forest types (Fig. 12b).The two differences are both large in the summer.For the difference between the GLOBCARBON LAI and the GLOBMAP LAI T , it ranges from 0.0 to −0.5 from November to April, while with values of up to −0.8 for ENFs, −1.3 for DNFs and DBFs in July and August.For the difference between the GLOBCARBON LAI and the GLOBMAP LAI o , it is smaller from November to April with a mean difference within 0.4, while it is larger for needleleaf forests during the summer, with mean differences up to approximately 1.0 for ENFs, 0.8 for DNFs and 0.4 for DBFs from June to August.
The forest understory vegetation and soil background vary with forest type, site and season.Notable uncertainties would be introduced if a constant background was used in the LAI retrieval.The GLOBCARBON LAI products remove partial effects of the forest background, leading to the differences between the GLOBCARBON LAI products and the GLOBMAP LAI T as well as the LAI o .These differences vary with the forest type, geographic location and season.Due to the different roles of the overstory and understory layers, uncertainties will be introduced into the estimation of the forest water and carbon cycle (Chen et al., 1999).This problem also exists for other LAI products which do not separate these two layers for forests.Pixel-specific seasonal forest background reflectivity helps to account for the effects of the forest background on the LAI retrieval.

Discussion
Many factors affect the quality of the derived forest overstory and understory LAI, such as LAI algorithm and inputs, mainly including MOD09A1 land surface reflectance, MISR forest background reflectance and clumping index.The MISR forest background reflectivity was used to estimate the forest understory LAI directly and correct for the effects of the background in the MODIS observations before the retrieval of the overstory LAI.Thus, the uncertainties in the MISR forest background reflectivity will be propagated into the results.It is pointed out that the tree architectural parameters used in the four-scale model and the land cover may affect the retrieved forest background reflectivity (Jiao et al., 2014), which will also affect the derived forest overstory and understory LAI.Additionally, the overstory and understory true LAI was converted from effective LAI retrievals with the global clumping index map derived from the MODIS BRDF product (He et al., 2012).The quality of this clumping index map would also affect our results.The reliability of the retrieved LAI is different for various forest types and seasons.Generally, it is difficult for satellite sensors to capture signals of the understory through the dense canopies.The reflectivities of the forest background tend to be unrealistic for dense canopies with LAI greater than 4, so the understory LAI is more reliable for sparse forest, and the retrievals in the spring and autumn should be more reliable than those in the summer.In addition, it is easier for optical remote sensing to capture signals from understory vegetation for coniferous forests than for broadleaf forests due to the more clumped foliage of the needleleaf forests.In fact, the structure of forest is very complicated.The overstory and understory may be composed of several sub-layers, and the young forests' canopy may be even continuous from ground to canopy top.In order to characterize the forest vertical structure from remote sensing data at a global scale, the forest is simplified to two layers, i.e., overstory and understory layers.Specifically, for EBF, its vegetation structure is much more complex than any other forest type, making it difficult to clearly separate the canopy into overstory and understory.In addition, in the forest background reflectivity algorithm, EBFs share tree architectural parameters (such as stand density, tree height and tree stem diameter) with DBFs, which may introduce significant uncertainties into the background reflectivity in the tropical zones due to large differences be-tween these two forest types.To avoid the influence of unconvincing forest background reflectivity, the EBFs were excluded in this study.
Forest understory has usually large species variation and heterogeneous spatial distribution.In addition, it is not entirely independent from the tree canopy since changes in canopy closure or tree layer LAI will lead to a change in the species composition and green LAI of ground vegetation (Rautiainen and Heiskanen, 2013).Generally, although the composition of understory is complex and site dependent, the typical species are shrubs, grasses and other herbaceous plants, mosses and lichens (e.g., Deering et al., 1999;Maeno and Hiura, 2000;Peltoniemi et al., 2005;Liang et al., 2012;Ryu et al., 2014;Qi et al., 2014;Nikopensius et al., 2015).In this paper, the understory LAI is estimated by averaging the retrievals based on the GLOBCARBON LAI algorithm for shrubs and grasses/crop/other non-forest vegetation.Thus, the complicated understory composition is simplified as shrubs and non-forest, non-shrub vegetation.The retrieved understory LAI may be affected by the uncertainties of the GLOBCARBON LAI algorithm for non-forest biomes.In addition, this simplification will also introduce uncertainties if the composition were not half shrubs and half other non-forest vegetation.Fortunately, the retrievals based on the algorithms for these two biomes are not quite differ- ent (Sect.3.1).Therefore, this simplification enables us to characterize the vertical structures for global forest from optical remote sensing observations.Moreover, lidar provides another powerful tool for monitoring forest vertical structures.For example, GLAS spaceborne waveform lidar data have been used to extract forest vertical foliage profiles over the United States at the footprint level (Tang et al., 2016).A combination of lidar and optical remote sensing may improve the forest structure monitoring.
It takes 9 days for MISR to acquire global coverage due to its relative narrow swath width of approximately 360 km.With influences of this low temporal resolution, occurrence of clouds, topographic obscuration and failed aerosol retrievals, a large proportion of missing data and invalid retrievals is included in the MISR land surface products (Liu et al., 2012b).As a result, there are large numbers of invalid retrievals in the MISR forest background reflectivity (Jiao et al., 2014).To generate spatially coherent maps, we have to compose the monthly forest understory LAI using the retrieved results during the 11-year period from 2000 to 2010.Thus, the changes in the forest understory among years and within months, such as fire and other disturbances, are not considered in this paper.Recently, the background reflectivity has been retrieved from MODIS BRDF products with a temporal resolution of 8 days over selected sites (Pisek et al., 2012(Pisek et al., , 2016)).Such BRDF products have high temporal resolutions and should make it possible to account for the interannual and seasonal variations in the forest understory vegetation in the future.

Conclusions
In this paper, the forest LAI is separated into overstory and understory layers over the global deciduous broadleaf and needleleaf forest areas based on monthly 1 km MISR forest background reflectivity and MODIS land surface reflectance.The global monthly 1 km forest understory LAI (GLOBMAP LAI u ) and 8-day 500 m forest overstory LAI (GLOBMAP LAI o ) were retrieved.
Forest overstory is the dominant component of forest LAI, and its spatial and seasonal patterns are close to those of the forest total LAI.The forest understory is mainly found in the boreal forest zones in the northern latitudes (50-70 • N), where 84 % of valid LAI u retrievals are found in July.Significant seasonal variations of the LAI u are present in this region, with LAI up to 2-3 from June to August.Needleleaf forests have more LAI u values greater than 1.5 than broadleaf forests, and their mean value (1.06 for ENFs and 1.04 for DNFs) is also larger than those for DBFs (0.96) due to their more clumped foliage.The mean proportion of understory LAI to total LAI is greater than 30 % (36 % for DBFs, 34 % for DNFs and 31 % for ENFs), indicating that the effects of forest understory cannot be ignored in canopy parameter estimation and forest carbon cycle modeling.The retrieved results show better consistency with the ground data for forests with lower LAI, while it is hard to separate the understory signal from the overstory for dense and closed canopies.
Although uncertainties still exist due to the many factors discussed above, this work would help us better understand the seasonal patterns of forest structure, distinguish different responses of forest overstory and understory to climate change, evaluate the ecosystem functions and improve the modeling of the forest carbon and water cycles.It is difficult to validate the forest understory and overstory LAI directly due to the lack of field measurements with separation of vertical layers, especially of understory LAI.Further studies of the structural characteristics of the forests and use of field lidar are necessary to evaluate the remote sensing overstory and understory LAI datasets.

Figure 1 .
Figure 1.General flowchart for the GLOBMAP forest overstory and understory LAI separation algorithm.

Figure 2 .
Figure 2. Frequency of differences of derived understory LAI based on the GLOBCARBON algorithm between shrub and grass.

Figure 3 .
Figure 3. Differences of derived understory LAI based on the GLOBCARBON algorithm between shrub and grass in (a) the Northern Hemisphere and (b) the Southern Hemisphere.

Figure 4 .
Figure 4. Global monthly GLOBMAP forest understory LAI maps by averaging valid retrievals during 2000-2010 for each month.Because notable uncertainties in the retrieved understory LAI are introduced for EBFs in the tropical zones due to the unreliability of the forest background reflectivities in this biome, the results from EBFs are excluded.

Figure 5 .
Figure 5. Global monthly GLOBMAP forest overstory LAI maps in 2010, with DOY001 for January, DOY041 for February, DOY073 for March, DOY105 for April, DOY137 for May, DOY169 for June, DOY185 for July, DOY225 for August, DOY249 for September, DOY281 for October, DOY313 for November and DOY345 for December.The results from EBFs are excluded.

Figure 7 .
Figure 7. Time series of the GLOBMAP overstory (thick lines) and understory (thin lines) LAI for (a) DBFs, (b) ENFs and (c) DNFs in the Northern Hemisphere.

Figure 8 .
Figure 8.Comparison between the GLOBMAP LAI o and LAI u with the Kobayashi datasets (Kobayashi et al., 2010) over the larch forest region in eastern Siberia: (a) understory LAI and (b) overstory LAI.

Figure 9 .
Figure 9.Comparison of the GLOBMAP LAI u (LAI o ) and ground LAI data: (a) results for understory LAI at Prince Albert National Park, Canada (N = 7); (b) results for overstory LAI (N = 28).

Figure 10 .
Figure 10.Seasonal profiles of GLOBMAP LAI u (blue) values and their comparison with in situ understory NDVI (NDVI u , green) measurements and estimated understory LAI from field NDVI u (field LAI u , pink) in seven forest sites, including stands in PFRR, Sudbury, Hyytiälä herb rich, Hyytiälä xeric, Järvselja RAMI birch, Järvselja RAMI pine and Laegern.

Figure 11 .
Figure 11.Maps of the differences (GLOBCARBON minus GLOBMAP) between (a) the GLOBCARBON LAI and the GLOBMAP forest total LAI (sum of the overstory and understory LAI) and (b) the GLOBCARBON LAI and the GLOBMAP forest overstory LAI.The results for the EBFs are excluded.

Figure 12 .
Figure 12.Monthly mean differences between (a) the GLOB-CARBON LAI and the GLOBMAP LAI T (sum of the overstory and understory LAI) and (b) the GLOBCARBON LAI and the GLOBMAP LAI o over the global forest area in the Northern Hemisphere based on data from 2008 to 2010.

Table 1 .
Summary of field LAI measurements used for validation with the derived forest overstory and understory LAI.CI stands for clumping status.For biomes, DBF, DNF and BF stand for deciduous broadleaf forest, deciduous needleleaf forest and boreal forest, respectively.For CI, the value Y (N) means the clumping effects have (have not) been taken into account in the LAI measurement.

Table 2 .
Summary of field NDVI measurements used for comparison with the derived forest understory LAI.

Table 3 .
(Kobayashi et al., 2010)n of the GLOBMAP LAI o and LAI u with the Kobayashi datasets(Kobayashi et al., 2010)over the larch forest in eastern Siberia.
paper, and it was compared with the Kobayashi LAI u .The GLOBMAP LAI o is compared with the Kobayashi LAI o from 2008 to 2009.The GLOBMAP LAI u and LAI o were resampled to the spatial resolution of 1/112 • and converted to the same geographic coordinates as the Kobayashi dataset.For each pixel, the GLOBMAP LAI o and LAI u were matched with the Kobayashi datasets at the corresponding location and the nearest observational dates to perform a pixelto-pixel comparison.Only those pixels labeled as DNFs both in MCD12 and GLC2000 were used.Figure • in geographic coordinates.A new forest understory LAI was estimated from the forest background reflectivities on leaf appearance day, which was provided in Kobayashi datasets, from 2005 to 2009 using the algorithm in this