Comparison of EPIC-Simulated and MODIS-Derived Leaf Area Index (LAI) across Multiple Spatial Scales

Abstract

Modeled leaf area index (LAI) in conjunction with satellite-derived LAI data streams may be used to support various regional and local scale air quality models for retrospective and future meteorological assessments. The Environmental Policy Integrated Climate (EPIC) model holds promise for providing LAI within a dynamic range for input into climate and air quality models, improving on current LAI distribution assumptions typical within atmospheric modeling. To assess the potential use of EPIC LAI, we first evaluated the Moderate Resolution Imaging Spectroradiometer (MODIS) LAI product collections 5 and 6 (i.e., Mc5, Mc6) with in situ LAI estimates upscaled at four 1.0 km resolution research sites distributed over the Albemarle-Pamlico Basin in North Carolina and Virginia, USA. We then compared the EPIC modeled 12.0 km resolution LAI to aggregated MODIS LAI (Mc5, Mc6) over a 3 × 3 grid (or 36 km × 36 km) centered over the same four research sites. Upscaled in situ LAI comparison with MODIS LAI showed improvement with the newer collection where the Mc5 overestimate of +2.22 LAI was reduced to +0.97 LAI with the Mc6. On three of the four sites, the EPIC/MODIS LAI comparison at 12.0 km resolution grid showed similar weighted mean LAI differences (LAI 1.29–1.34), with both Mc5 and Mc6 exceeding EPIC LAI across most dates. For all four research sites, both MODIS collections showed a positive bias when compared to EPIC LAI, with Mc6 (LAI = 0.40) aligning closer to EPIC than the Mc5 (LAI = 0.61) counterpart. Despite modest differences between both MODIS collections and EPIC LAI, the overestimation trend suggests the potential for EPIC to be used for future meteorological alternative management applications on a regional or national scale.

1. Introduction

The magnitude and seasonal variation in leaf area play a significant role in determining atmospheric deposition of airborne pollutants and, as such, is a key variable within atmospheric deposition modeling [1]. Satellite remote sensing of leaf area index (LAI), with LAI defined as the one-sided green leaf area per unit ground area in broadleaf canopies and as one-half the total needle surface area per unit ground area in coniferous canopies [2], allows for forest canopy characterization over large areas at broad spatial scales. However, satellite-derived LAI products can be limited by obstructed atmospheric conditions yielding sub-optimal values, or complete non-returns. Ecological research has utilized satellite-derived LAI estimates for retrospective analysis, e.g., [3,4], however these data will never be available for prospective meteorological or alternative management applications. This study is a two-tiered effort to compare regional LAI estimates from the most recent algorithm change in the Moderate Resolution Imaging Spectroradiometer (MODIS) LAI product and to compare these satellite estimates to LAI generated from an agricultural model developed by the United States Department of Agriculture (USDA). Modeled LAI inputs hold promise in support of various regional and local scale air quality models for retrospective and future meteorological assessments.

LAI can be a critical determinant of regional air and water quality [1]. Leaf area impacts the carbon sequestration process from carbon emissions [5], the depositional removal of atmospheric pollutants [6], the creation of tropospheric ozone from biogenic emissions of volatile organic compounds [7], and the regulation of the change of global meteorological patterns over time [8]. Several process-based atmospheric models are highly sensitive to small deviations of LAI. One model that ingests LAI is the Multilayer Model (MLM), a local scale atmospheric deposition model that estimates ozone (O₃), water vapor, nitric acid (HNO₃), and sulfur dioxide (SO₂) fluxes across the Clean Air Status and Trends Network (CASTNET) [9,10,11]. The MLM inputs a generalized annual LAI time-series profile developed from periodic leaf-on and leaf-off field sampled optical point measurements of each plant species. This model was shown to be highly sensitive to LAI profile parameters extracted from these point measurements generating estimates greater than 60% for ozone and ±25% for sulfur dioxide and nitric acid [12]. The regional Community Multiscale Air Quality Model (CMAQ), developed in part by the United States Environmental Protection Agency (EPA) [8], is partially reliant on Pleim-Xiu LAI (PX LAI), an LAI model producing estimates based on: (1) deep soil temperature and (2) an LAI response function [13,14]. This LAI response function relies on land-cover (LC)-specific minimum/maximum LAI and soil temperature, resulting in known overestimation biases [15]. CMAQ, MLM, and other atmospheric and air-quality models are likely to improve with more detailed temporal and spatial LAI inputs.

LAI has been estimated both at regional and global scales using satellite remote sensing data from a variety of sensors including the Moderate Resolution Imaging Spectroradiometer (MODIS) on NASA’s Terra and Aqua satellite platforms. MODIS LAI data originally processed in 2000 (collection 3), has been reprocessed three additional times with each iteration (collections 4–6) improving LAI returns through improved inputs and refinements in the main algorithm [16,17]. The most recent collection 6 (release date: 2015) improved upon the previous collection 5 (release date: 2007) through an improved surface reflectance product and the replacement of a static LC product with new multi-year LC product, both inputs benefiting from improved spatial resolution (1.0 to 0.5 km) [16].

LAI has also been simulated by a number of land-surface models including the United States Department of Agriculture (USDA) Environmental Policy Integrated Climate (EPIC) model [18], the ORganizing Carbon and Hydrology In Dynamic EcosystEms (ORCHIDEE) model [19] and the Interactions between Soil, Biosphere and Atmosphere, CO2-reactive (ISBA-A-gs) model [20,21]. Comparing growing season hindcast LAI simulations and phenological patterns observed from satellite-derived LAI is useful in detecting inconsistencies in both satellite and modeled products, thus leading to improvement with both data streams. For example, length of growing season and leaf onset were compared between modeled and satellite-derived LAI in the Euro-Mediterranean region [22]. Gross primary productivity (GPP), LAI, and vertical optical depth (VOD) were simulated by a conceptual ecohydrologic model (EcoHydr) then assessed in terms of its ability to reproduce the observed satellite derived products of GPP, LAI, and VOD from MODIS LAI products [23].

Whether LAI is simulated or generated from satellite observations, in situ and algorithm version comparisons are invaluable in the assessment of these products. MODIS LAI has been validated with in situ measurements at regional and global scales since its inception in 2001 [24]. At a regional level, the evaluation of LAI generated from two land surface schemes in southwest China found one model overestimating LAI by a factor of two with another model displaying close agreement with in situ LAI [25]. Comparisons across multiple iterations of LAI algorithm improvements allow the modeling community to understand the effect of LAI perturbations on backcasting and forecasting outcomes. Two globally based studies compared collections 5 and 6 to track overall global mean LAI trends [26] and to validate improvements in LAI estimation across each of the MODIS LC types [16]. On a regional level, MODIS collection comparisons of collections 3–5 were made to track LAI differences in corn and soybean crop types in Illinois [27].

Comparison of EPIC LAI to MODIS-derived LAI and in situ LAI will allow for its assessment as a potential variable input into air quality models (i.e., CMAQ). The effect of erosion on agricultural productivity initially drove the development of the USDA EPIC model in the early 1980s [28]. The model allows simulation of many processes important in agricultural management, but more recent studies have expanded to consider both natural and managed shrubland and forests [29,30,31]. Simulation using a daily time-step can be performed over 100 years or more to represent relatively long-lived natural or managed forest stands [32]. Examples of EPIC use within forested systems include: (1) the assessment of light interception in forecasting water usage in Kenyan semi-arid hedge intercrop, monocrop, and hedge monoculture (Senna spectabilis cv. Embu) systems [33]; (2) the examination of forest management practices on the loading of nutrients and sediments and streamflow in nine small Texas watersheds [34] (http://blackland.tamu.edu/models/apex); (3) the estimation of nitrogen, phosphorous, and herbicide losses across three undisturbed control watersheds, three conventional clearcut watersheds, and three intensive clearcut watersheds [35], and most recently; (4) our study that investigates EPIC model performance for predicting forest stand canopy height and LAI across four forest complexes in Virginia and North Carolina [31].

The objective of this study is to compare year 2002 MODIS LAI (collections 5 and 6—i.e., Mc5, Mc6) with (1) in situ LAI estimates upscaled to four 1.0 km resolution research sites distributed over the Albemarle-Pamlico Basin (APB) in North Carolina and Virginia, USA; and (2) EPIC modeled LAI over a 3 × 3 grid of nine 12.0 km grid cells centered over each of these four research sites. This two-fold effort will provide regional MODIS LAI validation across the prior (Mc5) and most recent (Mc6) algorithm updates and investigate the plausibility for modeled leaf area index (LAI) inputs to be used in conjunction with satellite-derived LAI data streams to support various regional and local scale air quality models for retrospective and future meteorological assessments.

2. Methods

• Comparison 1: In situ LAI to MODIS LAI (Mc5, Mc6) at 1.0 km.

Purpose: To explore statistical relationship between upscaled in situ LAI measurements and the latest MODIS algorithm update (Mc6) and check improvements over the previous iteration (Mc5).

• Comparison 2: MODIS LAI (Mc5, Mc6) to EPIC LAI at 12.0 km.

Purpose: To compare MODIS (Mc5, Mc6) LAI aggregated to 12.0 km to EPIC LAI over a 3 × 3 grid centered over the four sites to examine compatibility for use of EPIC at modeling air quality futures.

2.1. Site Descriptions

This research spans two physiographic provinces (coastal plain and piedmont) within the APB located in central-to-northern North Carolina and southern Virginia (Figure 1,

Table 1. Location, region, and tree species distribution over four 1 km LAI research sites.

Site	Lat/Long	Elevation (m)	Region	Tree Species
Appomattox	37.22, −78.88	165–215	Upper Piedmont	HDWD (21.8%), PINE (67.6%), OV (10.6%)
Hertford	36.38, −77.00	8–10	Coastal Plain	HDWD (17.8%), PINE (72.3%), OV (9.9%)
Fairystone	36.77, −80.09	381–473	Upper Piedmont	HDWD (100.0%)
Umstead	35.86, −78.73	75–120	Lower Piedmont	HDWD (55.1%), PINE (44.9%)

Note: HDWD—“hardwood”, OV—“other vegetation”.

). Four LAI research sites within the APB were chosen for sampling based on forest composition and access (

Table 2. Forest composition and biometrics for the four LAI research sites.

	Appomattox	Hertford	Fairystone	Umstead
1 Species (Dominant Canopy)	131	131	832, 802, 833	802, 837, 833, 621
Species (Understory)	316, 800	316, 591, 611	316, 407, 693	800, 403
Age (years)	23	19	80	80
2 DBH (cm)	21.6	19	24.8	43.5
Height (m)	15.9	14.3	18.9	24.4
3 TPH	1250	1740	288.8	101.5

¹ Tree species codes from the Forest Health Monitoring Program, United States Forest Service: 131—loblolly pine, 316—red maple, 403—pignut hickory, 407—shagbark hickory, 591—American holly, 611—sweetgum, 621—yellow-poplar, 693—blackgum, 800—oak spp., 802—white oak, 832—chestnut oak, 833—northern red oak. ² DBH—diameter at breast height, ³ TPH—trees per hectare.

). These sites included: 1. EPA_Hertford (Coastal Plain)—primarily loblolly pine; 2. EPA_Fairystone (Upper Piedmont-Mountain)—primarily upland oak; 3. EPA_Umstead (Lower Piedmont)—primarily upland oak/yellow poplar; 4. EPA_Appomattox (Upper Piedmont)—primarily loblolly pine (Table 1 and Table 2, Figure 1). Two sites are managed timberland (EPA_Appomattox and EPA_Hertford), the other two sites are managed state lands located in Virginia and North Carolina (Note: each research site name hereafter is truncated as: Hertford, Fairystone, Umstead, Appomattox). Selection of these LAI research sites corresponded to the 1.0 km grid associated with the 1.0 km MODIS (MOD15A2) 8-day composite LAI (collection 5) product produced by the National Aeronautics Space Administration (NASA). A detailed description of the forest composition of these four sites is also found in an initial study comparing EPIC-modeled LAI with in situ LAI (Supplementary Materials S1–S8, Table S1) [31].

2.2. LAI Data (In Situ LAI and Upscaling, MODIS LAI, EPIC LAI)

2.2.1. In Situ LAI Measurements and Upscaling

Sampling data was collected within the quadrat, a 100 × 100 m grid (Figure 2A), delineated by five 100 m east-to-west (E-W) oriented transects separated 20 m apart in a north-south (N-S) direction. Sampling point locations (25) were interspersed between these transects, separated 20 m apart. Secondary sampling units (i.e., sub-plots) comprised of 2–3 intersecting transects 50 and/or 100 m in length which were distributed randomly within the LAI research site extent (Figure 2B) [31]. One quadrat was placed within the Appomattox, Hertford, and Umstead LAI research sites. Four quadrats were located within the Fairystone site to account for the variability of terrain and aspect. Quadrat locations were chosen randomly within four 0.5 km sections (NW, NE, SE, SW) within the 1.0 km research sites (Figure 3).

At each LAI research site, overstory and understory forest measurements (i.e., height, diameter, stocking, crown closure, LAI) were made using both fixed area and point sampling methods as described in [31]. LAI was estimated for both deciduous and coniferous forest canopies from two indirect optical measurements: (1) the Tracing Radiation and Architecture of Canopies analyzer (TRAC) optical sensor and (2) digital hemispherical photography (DHP) [31,36,37,38]. The modified Beer–Lambert light extinction function was used to solve for LAI [1,31].

TRAC and DHP measurements were recorded across all four research sites for the following dates: Fairystone (1 May 2002, 25 June 2002, 4 September 2002), Appomattox (5 March 2002, 23 May 2002), Umstead (23 April 2002, 21 October 2002), and Hertford (5 March 2002). Measurements were not taken in the thinned pine at Appomattox (5 March) and the unthinned pine at Umstead (9 May). Pine LAI for both dates was estimated by comparing relative differences of SR LAI between 2002 data collects at each site, then reducing/increasing a known in situ LAI value by this percent difference.

In situ LAI was upscaled to a 1.0 km resolution coincident with the MODIS grids to create LAI reference maps (LAI-RM) at all four LAI research sites for MODIS LAI comparisons (Mc5 and Mc6) (Figure 4). USGS three-band natural color orthoimagery and cloud-free 30 m Landsat 5 (L5) Thematic Mapper and Landsat 7 (L7) Enhanced Thematic Mapper+ data were acquired from the USGS EarthExplorer repository for all four research sites with selections as close to the in situ collection dates as possible. L5 and L7 satellites have been active since 1983 (L5) and 1999 (L7), with L5 decommissioned in 2013. Both satellites acquire spectral data (0.45–2.35 µm) at a 30 m spatial resolution over six bands, along with one thermal band (10.40–12.50 µm) at 120 m (L5) and 60 m (L7). L5 and L7 imagery was clipped to the 1.0 km research sites coincident with the MODIS LAI Mc5 and Mc6 grids. LC percentages were calculated via heads-up digitization of the natural color leaf-off high-resolution (0.61 m) USGS orthoimagery using ArcMap. The normalized difference vegetation index (NDVI) was calculated from available L5 or L7 data for the hardwood forest stand type and regressed against the in situ LAI values. If correlations were acceptable (r² > 0.30), then LAI would be predicted via the regression model, otherwise the in situ LAI would be used. For loblolly pine, we employed the simple ratio (SR) LAI algorithm [39], with the SR LAI values calibrated to the in situ LAI measurements, then upscaled to the 1.0 km resolution. This algorithm ingests either a L5 TM or L7 ETM+ derivative ban—the vegetation index, SR, to predict LAI. Managed P. taeda LAI was predicted with approximately 14% error when comparing SR LAI with field-measured LAI [39]. Finally, in situ LAI mean and standard deviation was weighted by LC class over the 1.0 km resolution to arrive at an overall mean and standard deviation for the LAI-RM. LAI from young forest regeneration estimated solely using the SR LAI algorithm due to a more dominant presence of loblolly pine (Appomattox site only).

2.2.2. MODIS LAI

Both the 1.0 km MOD15A2 (Mc5) and the 0.5 km MOD15A2H (Mc6) LAI products are composited over an 8-day period and both are derived from a three-dimensional radiative transfer model driven by an atmosphere corrected surface reflectance product, a LC product, and ancillary information on surface characteristics. The half-kilometer resolution of the Mc6 collection reduces the uncertainty associated with heterogeneous landscapes, where one LC (i.e., biome type) is assigned to each pixel [16]. MODIS daily surface reflectances, previously produced at 1.0 km, were replaced by a 0.5 km product, which improved product output through better aerosol retrieval and algorithm corrections. In addition, the LAI product was improved through replacement of the static LC input with new multiyear LC maps. Both Mc5 and Mc6 MODIS LAI output is determined via two separate pathways: the main and the back-up algorithm. The main and back-up algorithms ingest an eight biome-specific classification map at 1.0 km (Mc5) and 0.5 km (Mc6), an atmospherically corrected MODIS spectral reflectance product (1.0 km—Mc5; 0.5 km—Mc6), and ancillary data containing radiative transfer coefficient lookup tables, backup algorithm lookup tables and output variable properties. The eight biome types used to generate Mc5 and Mc6 LAI data include: (1) grasses and cereal crops; (2) shrubs; (3) broadleaf crops; (4) savannas; (5) evergreen broadleaf forests; (6) deciduous broadleaf forests; (7) evergreen needle leaf forests; and (8) deciduous needle leaf forests. Biome-specific lookup tables containing the most probable LAI values are developed from iterative runs of a three-dimensional canopy radiative transfer model. In the case of the main algorithm, possible solutions (i.e., retrievals) are cases where the differences between the radiative transfer modeled and MODIS observed reflectances are within the uncertainty of the observed reflectances [40]. Dimensionless uncertainty assigned to the red and near infrared MODIS bands are estimated at 20% and 5% (non-forests) and 30% and 15% (forests), respectively [41,42]. Thus, a 1.0 km MODIS cell is assigned the mean of the retrieved LAI distribution.

The MODIS LAI (MOD15A2H), collection (6) Level 4 data was downloaded from the NASA EarthData website (https://earthdata.nasa.gov) for the 2002 year. The native projection for this data is the Sinusoidal Projection. MOD15A2H data, an 8-day composite data set with 0.5 km pixel size, is delivered in 1200 × 1200 km² tiles in the Hierarchical Data Format. The algorithm chooses the “best” pixel available from all the acquisitions of the MODIS Terra sensor from within the 8-day period. MODIS LAI is processed over eight biome types identified by the MODIS LC product (MCD12Q1), where each 0.5 km pixel is classified according to its status as a land or non-land pixel. Non-land pixels include: perennial salt or fresh water, barren or sparse vegetation, perennial snow or ice, permanent wetlands/inundated marshlands, or urban/built-up. The entire APB was contained within tile: H-11, V-5. All 46 MODIS LAI scenes for 2002 were stacked then reprojected into the Universal Transverse Mercator projection zones 17 and 18 (WGS 84 datum). This process was repeated for 1.0 km MODIS LAI (MOD15A2) collection (5) data downloaded from the USGS Earth Resources Observation and Science Center.

MODIS LAI (Mc5 and Mc6) in its raw form is not functional for analysis unless the data are smoothed to compensate for missing data and bad retrievals. Due to snow, ice, or cloud issues, LAI products often have pixels with anomalous high (spikes) and low (drops) values. Thus, the Adaptive Savitzky–Golay (ASG) filter was used to smooth LAI time series for the APB. This filter is a weighted moving average based on simplified least-squares fit convolution. The filter optimizes the values of the length of the temporal window and the order of the polynomial to get the best match between observations and reconstructed values. This technique tends to preserve higher LAI values. ASG takes advantage of cloud flag ancillary data where rogue LAI data points are replaced by linearly interpolated values using adjacent points. ASG is computationally simple and can be used to reconstruct high-quality LAI time-series by setting only two parameters: the half-width of the smoothing window and the degree of the smoothing polynomial [43]. We used TIMESAT package [44] to fit a Savitzky–Golay smooth function to the data. Pixels with low data quality flags were assigned with zero weight in the data smoothing process, thus they had no impact in the smoothing function.

2.2.3. EPIC LAI

In concept, the EPIC model is run to estimate LAI across multiple, competing species. This model uses information regarding species, species age, and species density to simulate a multi-age, multi-species “stand” which is then assumed to exist on all forested land in a study area [31]. For large area studies, the modeled stand characteristics usually derive from sampled stand characteristics. As such, it is unlikely that the modelled stand will equate exactly to any real world stand within the study area, but from a statistical perspective, it should be representative of the population of stands within the forestlands and could serve to constrain the range of possible stand-level LAI estimates. Stand characteristics were estimated for use in EPIC across the four EPA LAI research sites within nine cells of a 3 × 3 grid, consisting of nine 12.0 km grid cells (Figure 4). We used the US Forest Service Forest Inventory and Analysis (FIA) Evalidator Tool (http://apps.fs.fed.us/Evalidator/evalidator.jsp) to provide a consistent and reproducible estimate of species density and area information for an 88.3 km² area (5.3 km radius circle) surrounding the validation areas of interest. The FIA program samples forest land on public and private lands on average one plot for every 2428 hectares of forest land. In theory, this averages approximately 5.9 plots for every 12.0 km grid cell provided 100 percent is forested land. Stand age and diameter distribution by species were extracted for each of the nine 12.0 km grid cells per site. FIA data typically are aggregated into species groupings [31] that determine model growth parameters. In our prior study, EPIC LAI was calibrated thorough modification of distributed plant parameter values. Modified and/or added tree species plant parameter values (i.e., shade tolerance, maximum canopy height, species geographic distribution, temperature sensitivity, and seasonal growth pattern) were initially derived from the USDA Silvics manual [45]. These initial values were then calibrated using in situ observations of stand-level LAI and maximum canopy height while preserving reasonable physical and process driven relationships.

Within EPIC, tree species are assessed for species maturity, then cross-checked with diameter and age class distributions extracted from the FIA database. FIA Evalidator results typically return species groupings, i.e., sycamore/pecan/elm, sugarberry/hackberry/elm/green ash, etc. If more than one species shows comparable maturity dates, then multiple simulations (i.e., cohorts) are run with each species in that grouping contributing LAI in an additive manner. Stem densities are then partitioned between these species’ groupings, initiated (i.e., stems planted) at the year that brings maturity at year 2002, the year of our evaluation. All simulations were run from a simulated weather generator for the nearest National Oceanographic Atmospheric Administration’s (NOAA) National Centers for Environmental Information (NCEI) Cooperative Observer Program (COOP) locations (

Table 3. Cooperative Observer Program (COOP) sites and proximity to the four LAI research sites.

Site	COOP Site	Elevation Diff (m) COOP-Site	Distance (km)
Appomattox	51,011	77.4	15.8
Hertford	374,456	29.6	37.9
Fairystone	515,300	−239	21
Umstead	377,079	27.9	8.2

). Parameters extracted from these weather simulations include, maximum temperature (°C), minimum temperature (°C), and precipitation (mm). EPIC can be set to run both with climate generated variables or observed meteorological data, or both within a simulation.

EPIC requires stand establishment as a primary input to initiate stand development over the lifespan of the particular species. Other operational variables can be added including disturbance (i.e., ‘thinning’), harvest, or prescribed burn. Detailed multilayer soil profiles, commonly used simulation inputs for biogeochemical models, were not available for any of the four research sites. Instead we used a representative soil profile at all locations, the Cecil soil profile. Soil parameters for this profile resided within the Baumer database from the USDA, National Resources Conservation Service (NRCS) Soils Laboratory. The EPIC soil datasets were built to represent the sample point soils selected for the 1997 USDA Natural Resources Inventory (NRI) data points.

2.3. Comparisons: 1.0 km LAI-RM v. MODIS LAI and 12.0 km MODIS LAI v. EPIC LAI

2.3.1. Comparison 1: LAI-RM to MODIS LAI (Mc5; Mc6) at 1.0 km

We first compared LAI-RM (1.0 km resolution) to the corresponding MODIS LAI product (collection 5 (i.e., Mc5) [1.0 km resolution] and collection 6 (i.e., Mc6) [0.5 km resolution]). MODIS LAI Mc5 was compared one to one with the LAI-RM, whereas the Mc6 product was averaged over four 0.5 km cells within each 1.0 km area. MODIS variability was estimated temporally bracketing LAI values ± 1 MODIS LAI date (i.e., 8 days prior and post) on each side of the comparison date.

2.3.2. Comparison 2: MODIS LAI (Mc5; Mc6) and EPIC LAI at 12.0 km

To compare EPIC-derived LAI with MODIS LAI a 3 × 3 grid of nine 12.0 km resolution grid cells were created in ArcMap (i.e., ‘Fishnet’) centered about the four LAI research sites. The 12.0 km resolution was chosen to correspond with current air quality grid resolution for CMAQ. Spatial aggregation yielded LAI zonal statistics (mean, maximum, minimum, range, standard deviation) across the nine 12.0 km grid resolution cells for the smoothed MODIS LAI (Mc5 [144 pixels] and Mc6 [576 pixels]). Data quality flags were also extracted for both MODIS collections for ASG filter smoothing. EPIC-derived LAI was estimated at this same 12.0 km resolution for comparison. Notched box whisker plots were created to compare the medians of each data type at each date sampled. This informal graphical test displays a ‘notch’ where notch overlap signifies strong evidence (95% confidence) that the medians do not differ [39].

3. Results

3.1. Comparison 1: LAI-RM to MODIS LAI (Mc5; Mc6) at 1.0 km

3.1.1. LAI-RM (1.0 km) Creation from In Situ LAI

TRAC and DHP measurements were recorded and LAI calculated for all dates across all four research sites, then compared to Landsat NDVI (hardwood) and SR LAI (pine). The dates of in situ collection and Landsat imagery were coincident within a few days to a week except for two dates, one at Hertford and the other at Fairystone (

Table 4. In situ LAI, Landsat NDVI, MODIS collection dates.

	In Situ	Landsat	MODIS
Appomattox	5 March 23 May	5 March 24 May	5 March 23 May
Hertford	5 March	18 February	5 March
Fairystone	1 May	29 April	30 April
Umstead	25 June 4 September 9 May 21 October	25 June 6 October 1 May 23 October	25 Jun 6 September 9 May 24 October

). All four research sites, across all dates, exhibited poor relationships (r² < 0.30) between NDVI and in situ LAI for hardwood, therefore in situ hardwood LAI was selected for upscaling to 1.0 km. A comparison SR LAI and in situ LAI in the pine LC (thinned and unthinned) showed no difference in the means except for Umstead (24 October—unthinned) and Appomattox (5 March—unthinned) sites. In both cases, in situ LAI rather than SR LAI was used for upscaling. The four research sites (

Table 5. 1.0 km upscaled in situ LAI across four research sites (HDWD—hardwood, OV—other vegetation).

			In Situ LAI			LAI-RM
Date	EPA Site		n	$\overline{\mathit{x}}$	σ	Area (%)	$\overline{\mathit{x}}$	σ
23 May	Appomattox	HDWD	3	2.87	0.82	21.8	0.63	0.18
		PINE Unthinned	21	1.88	0.20	49.6	0.93	0.10
		PINE Thinned	8	1.01	0.19	18	0.18	0.03
		*** OV	0	0.38	0.27	10.6	0.04	0.03
						Total	1.78	0.34
5 March	Appomattox	HDWD	NA	0.00	0.00	21.8	0.00	0.00
		PINE Unthinned	24	1.49	0.14	49.6	0.74	0.07
		PINE Thinned	0	** 0.81	0.14	18	0.15	0.02
		*** OV	0	0.00	0.00	10.6	0.00	0.00
						Total	0.88	0.10
5 March	Hertford	HDWD	NA	0.00	0.00	17.8	0.00	0.00
		PINE Unthinned	20	1.73	0.16	72.3	1.25	0.12
		AG	0	0.00	0.00	9.9	0.00	0.00
						Total	1.25	0.12
30 April	Fairystone	HDWD	54	1.99	0.34	100	1.99	0.34
25 June	Fairystone	HDWD	84	2.54	0.38	100	2.54	0.38
4 September	Fairystone	HDWD	59	2.67	0.46	100	2.67	0.46
9 May	Umstead	HDWD	21	2.32	0.49	55.1	1.28	0.27
		PINE Unthinned	32	* 3.41	0.75	44.9	1.53	0.34
						Total	2.81	0.61
23 October	Umstead	HDWD	25	2.81	0.58	55.1	1.55	0.32
		PINE Unthinned	32	3.69	0.06	44.9	1.66	0.03
						Total	3.21	0.35

Note: * LAI estimated by comparing SR LAI change between the 23 October and 9 May estimates (7.49% decrease). ** LAI estimated by comparing SR LAI change between the 23 May and 5 March estimates (19.8% decrease). *** OV—2–3 year pine regeneration—used SR LAI to estimate LAI.

) varied in LC homogeneity over the 1.0 km with Appomattox exhibiting the most diverse ‘within’ LC types (i.e., natural pine, unthinned planted pine, pine thinned every few rows) and ‘between’ LC types (i.e., pine, hardwood, other vegetation). Fairystone was at the other end of the spectrum, with one LC type, hardwood, accounting for 100% of the 1.0 km area (Table 5).

3.1.2. Comparison of MODIS Mc6 0.5 km and MODIS Mc5 1.0 km to LAI-RM 1.0 km

Comparing LAI-RM (1.0 km) with the corresponding 1.0 km Mc5 across five of the eight dates showed Mc5 LAI values 1.0 to 2.4 times larger than that of the in situ upscaled LAI values (

Table 6. MODIS LAI (Mc5 and Mc6) comparisons to LAI-RM.

Site	MODIS Date	LAI-RM	LAI Mc5	LAI Mc6	Mc5—LAI-RM	Mc6—LAI-RM
Hertford	5 March	1.25	1.13	1.63	−0.12	0.38
Appomattox	5 March	0.88	1.93	0.88	1.05	0.00
Fairystone	30 April	1.99	5.23	2.52	3.24	0.53
Umstead	9 May	2.81	1.93	4.55	−0.88	1.74
Appomattox	23 May	1.78	6.40	4.13	4.62	2.35
Fairystone	25 June	2.54	5.33	3.29	2.79	0.75
Fairystone	6 September	2.67	4.57	3.32	1.90	0.65
Umstead	24 October	3.21	6.40	1.83	3.19	−1.38
	Mean	2.14	4.12	2.77	2.22	0.97
	SD	0.81	2.13	1.28	1.50	0.78

and Figure 5 and Figure 6). Mean Mc6 LAI exhibited better fit with LAI-RM estimates across these same seven dates with Mc6 exceeding 1.0 LAI units on only one date. Across all eight dates, the mean difference LAI-RM and MODIS LAI was 2.22 (SD = 1.50) for Mc5 and 0.97 (SD = 0.78) for Mc6.

3.2. Comparison 2: 12.0 km—EPIC LAI v. MODIS (Mc5/Mc6) LAI

3.2.1. MODIS Mc5 and Mc6 LAI Aggregated to 12.0 km

Mc5 and Mc6 LAI was processed for nearly 100 percent of the 1.0 km and 0.5 km pixels respectively within each 12.0 km grid cell for three of the four LAI research sites (Appomattox, Hertford, Fairystone) (Figure 7). The more urban dominated Umstead site, located in Wake and Durham counties, had a higher percentage of non-processed MODIS 1.0 km pixels across the nine grid cells ($\acute{x}$ = 43.1%, σ = 28%) as compared to the Mc6 pixels ($\acute{x}$ = 21.0%, σ = 25.1%). In general, the main radiative transfer algorithm was employed across early and late season growth to arrive at an LAI solution across all four research sites for both Mc5 and Mc6 MODIS products.

The mean LAI of the nine 12.0 km grid cells per site were similar for both the Mc5 1.0 km (n = 144) MODIS LAI and Mc6 0.5 km (n = 576) MODIS LAI products with Mc5 LAI exhibiting slightly higher LAI values than the Mc6 counterpart for the Appomattox site (3.0–4.7% LAI) (

Table 7. Mean LAI differences across each 3 × 3 grid of 12.0 km grid cells per site per date.

Site	Date	$\mathbf{Mc}5\overline{\mathit{x}}\left(\mathit{s}\right)$	$\mathbf{Mc}6\overline{\mathit{x}}\left(\mathit{s}\right)$	Mc5–Mc6	% Difference (Mc5 to Mc6)
Umstead	23 Apr	3.38 (0.54)	1.95 (0.78)	1.43	−42.3
Umstead	24 Oct	1.52 (0.23)	0.98 (0.42)	0.54	−35.5
Appomattox	6 Mar	0.85 (0.09)	0.81 (0.03)	0.04	−4.7
Appomattox	25 May	4.08 (0.36)	3.95 (0.40)	0.13	−3.2
Appomattox	5 Aug	3.94 (0.23)	3.82 (0.55)	0.12	−3.0
Fairystone	1 May	3.94 (0.19)	3.43 (0.43)	0.51	−12.9
Fairystone	26 Jun	4.46 (0.39)	4.42 (0.35)	0.04	−0.9
Fairystone	12 July	4.38 (0.44)	4.39 (0.34)	−0.01	0.2
Fairystone	6 Sep	3.86 (0.40)	3.47 (0.48)	0.39	−10.1
Hertford	6 Mar	0.88 (0.29)	0.96 (0.16)	−0.08	9.1
Hertford	15 Apr	1.51 (0.36)	1.95 (0.23)	−0.44	29.1
Hertford	18 Jun	3.04 (0.79)	3.3 (0.33)	−0.26	8.6
Hertford	5 Aug	3.04 (0.79)	3.7 (0.34)	−0.66	21.7
	Mean ($\overline{x}$)	2.92	2.83	0.09	13.46
	SD (s)	1.31	1.36	-	13.67

). Umstead Mc5 exhibited in a larger difference when compared to Mc6 LAI across the two dates (35.5–42.3% LAI). Fairystone exhibited larger Mc5 LAI values than Mc6 over three of the four dates (−0.9–12.9%). Hertford Mc6 LAI trended larger than Mc5 LAI (8.6–29.1). A grid-by-grid analysis for all four research sites (n = 113 12.0 km grids) over all dates showed Mc5 LAI exceeding Mc6 LAI over 71.0% of the grid cells by a mean LAI of 0.21 (

Table 8. 12.0 km grid-based mean differences per LAI research site over all Mc5 and Mc6 data collections.

Site	LAI Dataset	No. of Dates	No. of Grids	$\mathbf{Mean}\left(\overline{\mathit{x}}\right)\mathbf{Difference}$
** Umstead	Mc5	2	14	0.99
	Mc6		0	N/A
Hertford	Mc5	4	17	0.52
	Mc6		19	0.54
Fairystone	Mc5	4	30	0.41
	Mc6		6	0.64
Appomattox	Mc5	3	19	0.33
	Mc6		8	0.42
Total	Mc5	13	80	0.51
	Mc6		33	0.53
	Overall		113	0.21 *

Note: * LAI = 0.21 bias in MODIS Mc5 LAI; ** For Umstead, two of the possible nine 12.0 km grid cells were not included in LAI calculations due to dominance of urban LC in these grid cells; No.—‘number’.

). Only Hertford showed a majority of 12.0 km grid cells where Mc6 LAI exceeded the older collection.

3.2.2. EPIC Modeled LAI 12.0 km

Calibrated LAI values from our prior study were used to model LAI values for all four research sites. For example, calibrated pignut hickory LAI was used in place of the generic “Mixed Upland Hardwood” FIA forest stand type for the Appomattox, Fairystone, and Hertford sites. Apple was used as a surrogate for sassafras and persimmon on the Fairystone site. The Umstead site was unique in that it has been relatively undisturbed from harvest operations over the last 100 years. Therefore, maturity dates of loblolly pine were adjusted from 55 years to 90 years. We have provided an example of EPIC-modeled LAI on the Hertford research site (see Appendix A, Figure A1).

3.2.3. MODIS Mc5 and Mc6 LAI Comparison to EPIC LAI

At the Fairystone site, the weighted LAI means for both MODIS collections exceeded EPIC LAI across the four dates, with one date exception for Mc6 (Mc5: 4/4 dates exceedance, Mc6: 3/4 dates exceedance) (Figure 8 and Figure 9,

Table 9. Mean EPIC and MODIS (Mc5 and Mc6) LAI comparisons at 12.0 km scale within 3 × 3 grid (n = nine 12 km grid cells per date).

Fairystone	$1\mathbf{May}\overline{\mathit{x}}\left(\mathit{s}\right)$	$26\mathbf{June}\overline{\mathit{x}}\left(\mathit{s}\right)$	$12\mathbf{July}\overline{\mathit{x}}\left(\mathit{s}\right)$	$6\mathbf{September}\overline{\mathit{x}}\left(\mathit{s}\right)$
Mc5	3.94 (0.19)	4.46 (0.39)	4.38 (0.44)	3.86 (0.40)
Mc6	3.43 (0.43)	4.42 (0.35)	4.39 (0.34)	3.47 (0.48)
EPIC	2.72 (1.11)	3.21 (1.26)	3.28 (1.24)	3.26 (1.29)
Mc5-EPIC	1.22	1.25	1.10	0.60
Mc6-EPIC	0.71	1.21	1.11	0.21
Hertford	$6\mathrm{March}\overline{\mathit{x}}$(s)	$15\mathrm{April}\overline{\mathit{x}}$(s)	$18\mathrm{June}\overline{\mathit{x}}$(s)	$5\mathrm{August}\overline{\mathit{x}}$(s)
Mc5	0.88 (0.29)	1.51 (0.36)	3.04 (0.79)	3.59 (0.64)
Mc6	0.96 (0.16)	1.95 (0.23)	3.30 (0.33)	3.70 (0.34)
EPIC	1.64 (0.44)	1.72 (0.41)	2.97 (0.50)	2.97 (0.48)
Mc5-EPIC	−0.76	−0.21	0.07	0.62
Mc6-EPIC	−0.68	0.23	0.33	0.73
Appomattox	$6\mathrm{March}\overline{\mathit{x}}$(s)	$25\mathrm{May}\overline{\mathit{x}}$(s)	$5\mathrm{August}\overline{\mathit{x}}$(s)
Mc5	1.01 (0.09)	4.09 (0.40)	3.76 (0.24)
Mc6	0.81 (0.03)	3.95 (0.40)	3.82 (0.55)
EPIC	1.45 (0.60)	2.27 (0.71)	2.31 (0.73)
Mc5-EPIC	−0.44	1.82	1.45
Mc6-EPIC	−0.64	1.68	1.51
Umstead	$23\mathrm{April}\overline{\mathit{x}}$(s)	$24\mathrm{October}\overline{\mathit{x}}$(s)
Mc5	3.38 (0.54)	1.52 (0.23)
Mc6	1.95 (0.78)	0.98 (0.42)
EPIC	2.17 (0.97)	2.67 (1.20)
Mc5-EPIC	1.21	−1.15
Mc6-EPIC	−0.22	−1.69

). The weighted means were highest in this 100% hardwood site, varying under 0.10 LAI for both Mc5 (LAI = 1.34) and Mc6 (LAI = 1.27) (

Table 10. 12.0 km grid-based mean differences per LAI site over all EPIC and Mc5 LAI data collections (13).

Site	LAI Dataset	No. of Dates	No. of Grid Cells	Mean Difference
** Umstead	Mc5	2	8	1.18
	EPIC		6	1.49
Hertford	Mc5	4	20	0.64
	EPIC		16	0.55
Fairystone	Mc5	4	26	1.68
	EPIC		10	0.46
Appomattox	Mc5	3	21	1.43
	EPIC		6	0.81
Total	Mc5	13	75	1.28
	EPIC		38	0.72
	Overall		113	0.61 *

* Note: LAI = 0.61 bias in MODIS Mc5 LAI; ** For Umstead, two of the possible nine 12.0 km grid cells were not included in LAI calculations due to dominance of urban LC in these grid cells; No.—number.

and

Table 11. 12.0 km grid-based mean differences per LAI site over all EPIC and Mc6 LAI data collections (13).

Site	LAI Dataset	No. of Dates	No. of Grid Cells	$\mathbf{Mean}\left(\overline{\mathit{x}}\right)\mathbf{Difference}$
** Umstead	Mc6	2	4	0.74
	EPIC		10	1.63
Hertford	Mc6	4	21	0.66
	EPIC		15	0.48
Fairystone	Mc6	4	24	1.56
	EPIC		12	0.69
Appomattox	Mc6	3	21	1.38
	EPIC		6	1.00
Total	Mc6	13	70	1.19
	EPIC		43	0.88
	Overall		113	0.40 *

* Note: LAI = 0.40 bias in MODIS Mc6 LAI; ** For Umstead, two of the possible nine 12.0 km grid cells were not included in LAI calculations due to dominance of urban LC in these grid cells; No.—number.

). The LAI medians were significantly different between both collections and EPIC, except for 6 September and for the Mc6 collection only from 1 May (Figure 8). This trend of small median agreement between EPIC LAI and the MODIS LAI collections was seen at the other three research sites, with only Mc5 showing no significant differences for half of the LAI dates at Hertford (Figure 8). The Umstead site, like the Fairystone site, showed similar weighted mean differences for both MODIS LAI collections (Mc5 LAI = 1.31 and Mc6 LAI = 1.38) as did Appomattox with differences of 1.29 LAI for both Mc5 and Mc6. Hertford weighted LAI mean differences were the smallest among the four sites (Mc5 LAI = 0.60 and Mc6 LAI = 0.59). In summary, three of the four research sites showed similar weighted mean LAI differences (LAI 1.29–1.34), with both Mc5 and Mc6 exceeding EPIC LAI across most dates. Only Umstead showed EPIC LAI exceeding both MODIS LAI collections across both dates. For all four research sites, both MODIS collections showed a positive bias when compared to EPIC LAI, with Mc6 (LAI = 0.40) aligning closer to EPIC than the Mc5 (LAI = 0.61) counterpart (Table 10 and Table 11).

4. Discussion

4.1. Mc5, Mc6, and LAI-RM Comparisons (1.0 km)

The improvement of MODIS LAI primary algorithm outputs did align with [16] findings where total uncertainty fell within 1.0 LAI units in comparison of LAI-RM to MODIS LAI. The Mc6 improvements reduced the Mc5 overestimate (+2.22 LAI) to less than 1.0 LAI units (+0.97 LAI). This was also seen when evaluating median differences where only a little more than half of the dates showed no significant differences across all four research sites. An interesting outlier is seen comparing the Mc6 difference with the LAI-RM on the Appomattox pine-dominated site (67.6% pine). This site showed no LAI difference between both the Mc6 and LAI-RM estimates for March 5, while the late spring (23 May) LAI estimate was 2.35 larger for MODIS than for LAI-RM (Table 6). This exceedance is difficult to understand in relation to the in situ measurements of the differing LC types within the 1.0 km reference area. The pine in situ LAI more than doubled between the March and May dates due to added needle flushes and understory hardwood leaf-out. But even this increase in conjunction with the hardwood LC component, which represented 21% of the LC on the 1.0 km site yielding a LAI of 2.87, could not account for the Mc6 recorded 4.13 LAI.

Beyond error associated within the primary algorithm, two potential entry sources for error within the post-processing of MODIS LAI include the geometric reprojection from the native MODIS sinusoidal projection to the Landsat Universal Transverse Mercator projection and the temporal smoothing of the raw MODIS data. However, error introduced with the reprojection of the MODIS 1.0 and 0.5 km LAI data to the native Landsat projection has been reported as negligible [47,48]. Additionally, the ASG smoothing filter does add some error but is rated as one of the better linear interpolation methods in removing frequency noise [49].

Error in the creation of the LAI-RM may be introduced in the sampling strategy employed, in the optical measurements used within the modified Beer–Lambert light extinction function, in the analyst-effect in LC delineation [50], and in the correspondence of vegetation indices to ground measured LAI. In a prior study, we investigated variables ingested into the modified Beer–Lambert light extinction function. The largest sources of error identified were plant area index (PAI) and the woody-to-total ratio (W:T), variables subject to the thresholding method delineating sky from non-sky (PAI), and green vegetation from woody vegetation (W:T) [36]. It was understandable that we observed weak relationships between in situ LAI and vegetation indices derived from the red and near infrared bands of L5 and L7 imagery. The literature documents the NDVI:LAI correlation as suspect in forest stands exceeding a 3.0 LAI. Beyond LAI of 3.0, NDVI saturates i.e., an asymptotic increase with increasing LAI [51]. Saturation is a function of chemical and structural differences at the leaf and canopy levels (i.e., leaf orientation, pigmentation, woody-to-total ratio, canopy structure [tree height heterogeneity, crown classes present]) [48]. Because we did not have access to a good correlation in any of the hardwood stands, we also assumed our surrogate of the in situ LAI measured in a few areas was representative of the entire forest stand—an assumption that does not account for spatial heterogeneity of the forest canopy.

4.2. Mc5, Mc6, and EPIC Comparisons (12.0 km)

At the 12.0 km resolution grid cell scale Mc5 showed larger LAI values than Mc6 across 8 of the 13 dates for the four research sites, with values differing by 13.5 percent over all research sites and dates. We reason these differences at the grid cell level were related to input changes within the new Mc6 collection. One significant change between both collections was an improved biome map with Mc6 LAI reducing errors in biome classification. However, these differences may not have been recognized at the 12.0 km resolution grid cell level based on averaging of any propagated errors.

Across all four research sites and all 13 dates, EPIC and MODIS LAI (Mc5 and Mc6) showed similar values when evaluating directional bias, with MODIS exceeding EPIC LAI by 0.61 and 0.40 LAI, respectively. However, there were wide variations between both MODIS LAI collections and EPIC LAI on a site and date basis. This would preclude interchanging EPIC for MODIS on a site-specific (i.e., 36 × 36 km²) basis based on the lack of a directional bias. However, on a larger regional level, EPIC and MODIS might be interchangeable, appropriate for regional and national modeling of air quality. In a recent study, MODIS satellite products that included both LAI and the fraction of photosynthetically active radiation (FPAR) were integrated into the WRF-CMAQ modeling system and showed greater error and bias with temperature but reduced error and bias with the LAI and FPAR inputs [15].

The primary areas on uncertainty within the EPIC model include characterization of the vegetation component and the ingestion of the chosen climatological time series. The coarse distribution of FIA sampling points within each 3 × 3 grid (or nine 12.0 km cells) affected EPIC LAI outputs. This was seen by running multiple iterations using differing tree species associations when presented with an FIA tree not previously calibrated in the field. These FIA limited points (usually 2–3 per 36 km²) were assumed to completely describe all tree species associations throughout the entire 3 × 3 grid. Forest species groups may be predicted with greater accuracy using remotely sensed imagery in combination with other spatially continuous geospatial data. The USFS is developing a dataset that combines in situ FIA data with STATSGO (soil bulk density, depth to bedrock, soil permeability, soil porosity, available water capacity, soil plasticity, rock volume, soil types, soil texture, and soil pH), MODIS vegetation indices (Enhanced Vegetation Index [EVI], NDVI), prism temperature and precipitation data, and various other datasets to arrive at a distribution of 28 forest type groups across the continental United States. (https://databasin.org/datasets/7fa4ddf2a1c447469dbe7e98f81de814). This 250 m product may provide a useful input to parameterize EPIC LAI simulations at regional or national scales.

Another source of variability in the EPIC simulations that was highlighted during the generation of the stand simulations is the importance of the climatological time series employed. Daily weather conditions were “generated” using statistics derived from time series (e.g., 50 years or more) of observed weather conditions. Such long time-series are used to ensure that extreme events across a range of return periods e.g., multi-year droughts, heat waves, and floods are represented. Such extreme events can damage or destroy forest stands. EPIC simulations, however, assume these events occur everywhere within the study area and so occasionally a simulated stand will fail to survive, and so no viable stand will be simulated for the entire study area, a situation which is unrealistic. This finding suggests that future applications should either employ higher resolution (spatial and/or temporal) weather information or a spatial algorithm is needed to better represent the spatial extent of potentially devastating extreme events on long-lived stands.

5. Conclusions

The Mc6 LAI product was improved when comparing the differences from in situ LAI (LAI-RM) from those recorded by Mc5. LAI-RM comparison with MODIS LAI showed improvement with the newer collection where the Mc5 overestimate of +2.22 LAI was reduced to +0.97 LAI with the Mc6. The EPIC LAI simulations compared favorably to both MODIS LAI collections across all three of four research sites and all 13 dates when evaluating the overall mean difference. On three of the four sites, the EPIC/MODIS LAI comparison at 12.0 km resolution grid showed similar weighted mean LAI differences (LAI 1.29–1.34), with both Mc5 and Mc6 exceeding EPIC LAI across most dates. For all four research sites, both MODIS collections showed a positive bias when compared to EPIC LAI, with Mc6 (LAI = 0.40) aligning closer to EPIC than the Mc5 (LAI = 0.61) counterpart. However, there was larger variability and directionality when viewing these differences on a site by site basis. EPIC LAI potentially could augment MODIS LAI data feeds at the EPIC coarser resolution as more forest types are calibrated and validated and as other inputs are improved. EPIC shows promise for forecasting future meteorological alternative management applications, where the potential of providing a dynamic range of LAI values from EPIC may supersede current LAI inputs found in many atmospheric models.