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

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.


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 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.  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.

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). Four LAI research sites within the APB were chosen for sampling based on forest composition and access (Table 2). 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 (Tables 1 and 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].  Note: HDWD-"hardwood", OV-"other vegetation".    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). 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    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: . 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. 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 2 > 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).  [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).

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 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 2 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),

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  (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 2 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 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.

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 2 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). 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. ). 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].  (Table 4). All four research sites, across all dates, exhibited poor relationships (r 2 < 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) 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).            (Figures 8 and 9, Table 9). The weighted means were highest in this 100% hardwood site, varying under 0. 10

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).

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) (Figures 8 and 9, Table 9). The weighted means were highest in this 100% hardwood site, varying under 0. 10

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].  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. 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 2 ) 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 2 ) 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.

Conclusions
The  Table S1: Forest stand structural attributes for 4 sites.
Author Contributions: J.S.I. coordinated this research, establishing cooperation from private and public entities for the collection of the in situ data. He has analyzed the MODIS-EPIC relationship and has written the majority of this paper. Y.S. ran the smoothing algorithm on the MODIS time-series data. He also wrote the section on the smoothing algorithm. E.C. provided the EPIC LAI simulations for this evaluation. She wrote the EPIC methods and contributed to the results. A.P. assisted in the in situ data collection and the writing of the in situ LAI collection techniques. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.