Albedo Impacts of Changing Agricultural Practices in the United States through Space-Borne Analysis

: Using the collocated Moderate Resolution Imaging Spectroradiometer (MODIS)-derived Bidirectional Reﬂectance Distribution Function (BRDF) with the U.S. Department of Agriculture’s National Agricultural Statistics Service’s Cropland Data Layer (CDL), the daily albedo of homogenous agricultural ﬁelds was derived for 51 common United States ﬁeld crops by wavelength, sky-type, day of year, crop, and hardiness zone from 2015–2018. This study suggests that crop growth can result in changes in reﬂectivity up to a factor of 2 at most wavelengths and is unique per crop type in timing and range. Additionally, broadband impacts were studied and found to be less conspicuous than the individual wavelengths, but still signiﬁcant. The results were used to evaluate a common method of cropland albedo estimation, normalized di ﬀ erence vegetation index (NDVI) as a proxy for albedo, and this method was found to have some signiﬁcant limitations dependent on wavelength and date. Finally, a database of surface albedo variations as a function of growing period is constructed for common ﬁeld crops to the United States (as well as additional land-cover types). This database can be used to aid both satellite remote-sensing applications and long-term weather modeling e ﬀ orts by providing a method for parameter adjustments based on crop driven albedo changes, including changes in cropland composition related to commodity markets and other external factors.


Introduction
Surface albedo, which is the ratio of outgoing to incoming irradiance at a given wavelength, is a key component to a range of scientific applications, such as the lower boundary conditions for the remote sensing of surface and atmospheric properties [1], climate studies [2], and atmospheric modeling parameters [3]. Changes in surface albedo over cropland occur at weekly, and seasonal scales due to crop growth. Yet, while seasonal changes in surface characteristics are accounted for in some applications, these weekly changes in surface albedo are less frequently considered. Still, significant variation in both timing and range of surface albedo are observed for crops from initial planting to harvest stages, which can be noted simply by the significant variation of color hues displayed over croplands during maturation (e.g., Figure 1). To maintain soil health, suppress pests and disease, and thereby maximize long-term productivity and profits, farmers routinely use cyclical crop rotations at a given location such that the crop planted (along with pest and disease pressure) varies on an annual to semi-annual basis. As a result of historically routine crop rotation practices and new rotational practices related to market conditions and United States agricultural policy, surface characteristics of cropland derived from prior years may not be representative of current or future conditions. Additionally, in-season surface albedo variation caused by both crop rotation and crop growth throughout the growing season are a non-negligible issue due to the large spatial extent of agricultural lands. For example, in the continental United States in 2019, farmland covered a total of 895 million acres [4] (47.4% of total land area), with 308 million acres of that in harvestable crops [5] (16.3% of total land area).
Several approaches have been attempted to account for changes in surface albedo resulting from plant growth. For example, Houldcroft et al. [6] defined albedo parameters for five vegetation categories (and the corresponding bare soil) worldwide using Moderate Resolution Imaging Spectroradiometer (MODIS)-derived albedo, calculating albedo as a combination of plant and soil reflectivity. Hsu et al. [7] used the normalized difference vegetation index (NDVI) to account for cropland surface albedo changes for satellite aerosol retrievals, whereas Zhang et al. [8] produced time-series forecasts for albedo of a specific area including croplands through empirical mode decompositions and neural networks. Still, for most prior research applications, surface albedo values have either been assumed to be invariant throughout the study period or estimated through the use of a static look-up table based on generic cropland albedo changes.
In this study, through collocation of the U.S. Department of Agriculture's (USDA) National Agricultural Statistics Service's (NASS) Cropland Data Layer (CDL) with the MODIS bidirectional To maintain soil health, suppress pests and disease, and thereby maximize long-term productivity and profits, farmers routinely use cyclical crop rotations at a given location such that the crop planted (along with pest and disease pressure) varies on an annual to semi-annual basis. As a result of historically routine crop rotation practices and new rotational practices related to market conditions and United States agricultural policy, surface characteristics of cropland derived from prior years may not be representative of current or future conditions. Additionally, in-season surface albedo variation caused by both crop rotation and crop growth throughout the growing season are a non-negligible issue due to the large spatial extent of agricultural lands. For example, in the continental United States in 2019, farmland covered a total of 895 million acres [4] (47.4% of total land area), with 308 million acres of that in harvestable crops [5] (16.3% of total land area).
Several approaches have been attempted to account for changes in surface albedo resulting from plant growth. For example, Houldcroft et al. [6] defined albedo parameters for five vegetation categories (and the corresponding bare soil) worldwide using Moderate Resolution Imaging Spectroradiometer (MODIS)-derived albedo, calculating albedo as a combination of plant and soil reflectivity. Hsu et al. [7] used the normalized difference vegetation index (NDVI) to account for cropland surface albedo changes for satellite aerosol retrievals, whereas Zhang et al. [8] produced time-series forecasts for albedo of a specific area including croplands through empirical mode decompositions and neural networks. Still, for most prior research applications, surface albedo values have either been assumed to be invariant throughout the study period or estimated through the use of a static look-up table based on generic cropland albedo changes.
Remote Sens. 2020, 12, 2887 3 of 19 In this study, through collocation of the U.S. Department of Agriculture's (USDA) National Agricultural Statistics Service's (NASS) Cropland Data Layer (CDL) with the MODIS bidirectional reflectance distribution function data, the variations of surface albedo during the crop growing season are studied for 51 field crops, seven spectral channels (ranging from visible to shortwave infrared), and nine different plant hardiness zones across the United States. From this collocation, a database was created, including both the 51 crop types focused on in this study as well as 49 additional land-cover types, for applications that require accurate estimations of daily surface albedo changes over cropland. Using this database, this study tackles three overarching questions: (1) what are the impacts of crop variations in growth and reflectivity on spectral surface albedo, and are these impacts consistent on a per-wavelength, crop type, day of year, and hardiness zone basis?; (2) can cropland NDVI values be used as a proxy for the variations in surface albedo over cropland caused by plant growth?; and, (3) are the broadband albedo changes due to crop growth significant for climate applications?

Materials and Methods
The changes in spectral albedo as a function of crop growth were studied over the continental United States during a four-year period (2015-2018). Four datasets, including USDA-NASS's CDL, USDA's Plant Hardiness Zones database, and MODIS-derived albedo and reflectivity were collocated and used for the study. The study's region and timing were chosen due to the availability and consistency of the CDL data, which allows for identification and inclusion of a large number of homogeneous single-crop locations overlaying MODIS pixels.

Moderate Resolution Imaging Spectroradiometer (MODIS) Bidirectional Reflectance Distribution Function (BRDF) Albedo Product
The MODIS Bidirectional Reflectance Distribution Function (BRDF)/Albedo Parameter level 3, 500 m gridded, dataset (collection 6) is used for estimating surface albedo at seven wavelengths (0.47, 0.56, 0.65, 0.86, 1.24, 1.64, and 2.13 µm). This dataset calculates the effective black-and white-sky albedo parameters over a 16-day moving window [9]. Through black-and white-sky albedo, the surface albedo can be calculated for specific solar zenith angle and atmospheric scattering ratio, from Lucht, et al. [10] as follows: where α is the total albedo, SR di f is the ratio of indirect (diffuse) to direct sunlight reaching the surface, α w is the white-sky albedo, and α b is the black-sky albedo. This dataset further includes quality control flags for both the overall pixel quality and the snow coverage of any given pixel within the area during the 16-day moving window. Throughout this study, only snow/ice-free "good" or "best" quality pixels were utilized to focus on changes seen to the land surface due to crop growth and/or producer management. The MODIS BRDF dataset additionally provides a calculated broadband albedo for three spectral ranges, visible (300-700 nm), infrared (IR, 700-5000 nm), and shortwave (300-5000 nm), derived from the reflectance of the individual channels as shown in Equation (2) [10]: where the α r is the albedo of the spectral range, α i is the albedo as calculated for each independent channel, and w i is the weighting factor for each channel. w i was calculated through computational means to derive the relative impact of each channel to the overall albedo.

Cropland Data Layer (CDL)
To calculate the impact of crop type on surface albedo, precise information about the land cover homogeneity within each MODIS pixel is required. In this study the NASS CDL provided  locations and crop identifications for all cropland and other land uses within the United  States at a 30 m resolution. The CDL breaks down the land cover at a 30 m resolution by crop type  or non-crop land use through the application of a selection of independent satellite measurements,  including RESOURCESAT-1 and Landsat, to derive the reflectance of the shortwave bands for CDL pixel classification. Land use is then defined through a supervised classification technique that incorporates USDA Farm Service Agency (FSA) ground truth in combination with the remotely-sensed reflection values [11]. Note that uncertainties exist in the CDL data, which is discussed in detail in Section 3.3. The CDL was then collocated with the MODIS dataset as described in Section 2.5 to establish locations of homogeneous crop fields used in this study.

Plant Hardiness Zone
To account for changes in growing seasons and their impacts on the growth cycle of each crop, the selected locations are separated into specific hardiness zones (HZ) using the USDA's Plant Hardiness Zone (PHZM) defined zones [12]. The PHZM zones are split into half-zone divisions of 2.8 • C defined through the calculation of the 30-year average extreme minimum temperatures of 7983 North American weather stations. These are then interpolated, with elevation impacts incorporated, into the final map [13]. However, in order to streamline the categories in the final database, this study combines the PHZM half-zone divisions into their full zone (5.6 • C) counterparts.
By accounting for the growing season length and temperature, the calculations of crop albedo can be grouped by commonality of growth cycles. This is possible as the local climate largely defines both the planting and harvest dates, along with the growth rate of particular crops, resulting in common growth patterns within a given HZ. Grouping by HZ could potentially lead to reduction in uncertainty due to crops planted and harvested at similar times having a similar growth profile throughout the year.

MODIS Reflectance Data
To aid in the application of using NDVI as a proxy for per-channel albedo, level 2 MODIS surface reflectance (MOD09 + MYD09) [14] was utilized. Data were extracted from both Aqua and Terra satellites at a 500 m resolution for the 2015-2018 seasons and cloud-cleared using the internal cloud and cloud shadow masks. The NDVI calculation utilizes the red (645 nm center) and near infrared (860 nm center) channels to determine the vegetation index as follows from Rouse Jr, et. al. [15]: where REF NIR is the reflectance of the near infrared MODIS channel (860 nm) and REF Red is the reflectance of the red MODIS channel (650 nm). The resulting unitless index ranges from −1 to 1. NDVI allows for the analysis of plant growth and activity due to the increased absorption of the red channel and increased reflectance, relative to the underlying soil, of the near infrared channel of the chlorophyll active within plants. NDVI values >0.4 are generally associated with levels of plant growth with especially high values (>0.8) possible during times of peak growth of field crops. Finally, for the investigation into NDVI as a proxy of albedo, both the NDVI and albedo pixels were masked using the MODIS annual land cover in place of the CDL layer, allowing for both single and multi-crop MODIS pixels to be included in this aspect of the study. Note that as a potential source of error in this aspect, the MODIS satellite has a pointing uncertainty of around 50 m [16], which may lead to unforeseen issues when attempting to collocate over specific fields. Further potential uncertainties are discussed in Section 3.3.

Collocation of CDL and MODIS Data
Calculating the seasonal change in albedo of specific crops in this study requires collocation of the CDL with the MODIS gridded BRDF dataset. Figure 2 shows the flow chart of the collocation steps. Colocation was performed by locating the center of each BRDF grid pixel and locating the corresponding CDL pixel found at that location. Due to the resolution difference, (30 m CDL vs. 500 m MODIS) an appropriate buffer area around the BRDF pixels was selected extending the CDL pixel analysis zone beyond the MODIS pixel. All CDL pixels within the area + buffer were checked for single crop homogeneity, with the overlaying BRDF pixel discarded from the study if any pixel within the area + buffer was classified as a different crop or land cover from the central CDL pixel. As a result, all MODIS BRDF pixels used in this study are considered to be comprised entirely of a single crop as interpreted by the CDL. For the remaining, single crop, BRDF pixels were grouped by crop and HZ, and the mean albedo, standard deviation of that albedo, and raw pixel counts for each day of the year were calculated from the merged four-year BRDF dataset.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 18 single crop homogeneity, with the overlaying BRDF pixel discarded from the study if any pixel within the area + buffer was classified as a different crop or land cover from the central CDL pixel. As a result, all MODIS BRDF pixels used in this study are considered to be comprised entirely of a single crop as interpreted by the CDL. For the remaining, single crop, BRDF pixels were grouped by crop and HZ, and the mean albedo, standard deviation of that albedo, and raw pixel counts for each day of the year were calculated from the merged four-year BRDF dataset. As the result of colocation, 20,345,626 MODIS pixels across all land-cover types were found that matched the criteria over the four year period (roughly 5,000,000 per year); which are shown below grouped by land cover type (Figure 3a), by HZ (Figure 3b), and the selected points of several common US crops ( Figure 4) for reference. Using these reference points both mean and standard deviation of albedo for each crop, HZ, day of year, wavelength (including broadband), and sky-type combination were calculated. Additionally, means and standard deviations of the albedo for the remaining homogeneous, but non-cropland covers were derived and included in the final database, but not focused upon in this study. As the result of colocation, 20,345,626 MODIS pixels across all land-cover types were found that matched the criteria over the four year period (roughly 5,000,000 per year); which are shown below grouped by land cover type (Figure 3a), by HZ (Figure 3b), and the selected points of several common US crops ( Figure 4) for reference. Using these reference points both mean and standard deviation of albedo for each crop, HZ, day of year, wavelength (including broadband), and sky-type combination were calculated. Additionally, means and standard deviations of the albedo for the remaining homogeneous, but non-cropland covers were derived and included in the final database, but not focused upon in this study. 6

Results
The black-and white-sky albedo are studied as functions of crop type and growing period over the four-year study period in this section. Again, the goal of this effort is to construct a database providing for the fast derivation of albedo changes over cropland through the growing period, examples of the seasonal variations are shown in Table 1. The database includes the mean and standard deviation of the albedo, at each MODIS visible and near IR wavelength, including broadband, for both sky types, for every day of the year for a total of 51 field crops plus 49 non-fieldcrop land covers. Also included in the dataset are ancillary parameters such as valid MODIS pixel counts. The dataset is included as a supplement to this paper. Examples of the albedo changes as a function of growing days for visible, near IR, and shortwave infrared (SWIR) channels, as well as the broadband SW spectrum, are shown in the following sections for maize, soybean, spring wheat and cotton. These are four of the most commonly planted spring to fall cycle crops in the US, and are shown as examples; (Maize with 37.1 mil ha, soybeans with 32.4 mil ha, spring wheat with 5.0 mil ha, and cotton with 5.5 mil ha [5]).

Results
The black-and white-sky albedo are studied as functions of crop type and growing period over the four-year study period in this section. Again, the goal of this effort is to construct a database providing for the fast derivation of albedo changes over cropland through the growing period, examples of the seasonal variations are shown in Table 1. The database includes the mean and standard deviation of the albedo, at each MODIS visible and near IR wavelength, including broadband, for both sky types, for every day of the year for a total of 51 field crops plus 49 non-field-crop land covers. Also included in the dataset are ancillary parameters such as valid MODIS pixel counts. The dataset is included as a supplement to this paper. Examples of the albedo changes as a function of growing days for visible, near IR, and shortwave infrared (SWIR) channels, as well as the broadband SW spectrum, are shown in the following sections for maize, soybean, spring wheat and cotton. These are four of the most commonly planted spring to fall cycle crops in the US, and are shown as examples; (Maize with 37.1 mil ha, soybeans with 32.4 mil ha, spring wheat with 5.0 mil ha, and cotton with 5.5 mil ha [5]).     Figure 5 shows the change in black-sky surface albedo as a function of the growing period for maize, soybeans, spring wheat, each in HZ4, and cotton, in HZ8, for the three spectral MODIS channels (470 nm, 555 nm, 645 nm). The solid lines show the daily mean albedo for the seven wavelengths, while the shaded area shows the standard deviation. For illustration purposes, for maize, soybeans, and spring wheat, HZ4 is shown as the spatial overlap of these three crops is high in this HZ, allowing for direct comparison. However, as displayed in Figure 4, cotton does not appear in HZ4 and as such cotton in HZ8 is included for reference but does not spatially overlap with the other three crops.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 18 Figure 5 shows the change in black-sky surface albedo as a function of the growing period for maize, soybeans, spring wheat, each in HZ4, and cotton, in HZ8, for the three spectral MODIS channels (470 nm, 555 nm, 645 nm). The solid lines show the daily mean albedo for the seven wavelengths, while the shaded area shows the standard deviation. For illustration purposes, for maize, soybeans, and spring wheat, HZ4 is shown as the spatial overlap of these three crops is high in this HZ, allowing for direct comparison. However, as displayed in Figure 4, cotton does not appear in HZ4 and as such cotton in HZ8 is included for reference but does not spatially overlap with the other three crops. As shown in Figure 5, significant changes in black-sky cropland albedos in the visible spectrum, particularly with respect to the 645 nm channel, are clearly observable. For maize fields, black-sky albedo is reduced from 0.1 at the beginning of the planting season (day 100) to a low of 0.03 during the peak growing season (day 200), representing a 70% relative decrease in black-sky albedo possibly due to plant photosynthesis. Following the peak growth period, the black-sky albedo for maize fields increases to 0.15 to 0.2 during and after the harvesting period (after day 275), as plant photosynthesis activity diminishes during crop maturity and the crop is ultimately removed from the surface.

Variations in Cropland Albedo at Visible, Near Infrared (NIR), and Shortwave Infrared (SWIR) Channels Due to the Crop Growth Cycle
Near identical patterns in black-sky albedo are found for soybean fields compared with maize fields, suggesting that the two crops share similar growth and photosynthesis activity cycles. In comparison, spring wheat fields are planted earlier and mature at a much faster rate relative to maize and soybeans. Consequently, the minimum black-sky albedo is centered around late June (day 175), which is around 25 days earlier than soybean or maize fields (day 200), and the change is sustained for a much shorter period of time.
The black-sky albedo change for cotton is shown in Figure 5 for HZ8, which again is chosen as a representative zone as cotton is not commonly found in HZ4. Note that cotton fields in the United States are planted around mid-May (day 150) and harvested in late October (day 290). Greater variance and a less significant overall change was found during the peak growing season, with the black-sky albedo of the 645 nm channel changing from 0.15 to 0.10. This lack of definable response could be a result of the agricultural practices around cotton; as cotton has a wider range of planting and harvesting dates than maize and soybean, as well as the extended growth cycle relative to the other US annual crops based on the USDA-NASS annual report [17]. Additionally, the physiology of As shown in Figure 5, significant changes in black-sky cropland albedos in the visible spectrum, particularly with respect to the 645 nm channel, are clearly observable. For maize fields, black-sky albedo is reduced from 0.1 at the beginning of the planting season (day 100) to a low of 0.03 during the peak growing season (day 200), representing a 70% relative decrease in black-sky albedo possibly due to plant photosynthesis. Following the peak growth period, the black-sky albedo for maize fields increases to 0.15 to 0.2 during and after the harvesting period (after day 275), as plant photosynthesis activity diminishes during crop maturity and the crop is ultimately removed from the surface.
Near identical patterns in black-sky albedo are found for soybean fields compared with maize fields, suggesting that the two crops share similar growth and photosynthesis activity cycles. In comparison, spring wheat fields are planted earlier and mature at a much faster rate relative to maize and soybeans. Consequently, the minimum black-sky albedo is centered around late June (day 175), which is around 25 days earlier than soybean or maize fields (day 200), and the change is sustained for a much shorter period of time.
The black-sky albedo change for cotton is shown in Figure 5 for HZ8, which again is chosen as a representative zone as cotton is not commonly found in HZ4. Note that cotton fields in the United States are planted around mid-May (day 150) and harvested in late October (day 290). Greater variance and a less significant overall change was found during the peak growing season, with the black-sky albedo of the 645 nm channel changing from 0.15 to 0.10. This lack of definable response could be a result of the agricultural practices around cotton; as cotton has a wider range of planting and harvesting dates than maize and soybean, as well as the extended growth cycle relative to the other US annual crops based on the USDA-NASS annual report [17]. Additionally, the physiology of the cotton plant, which includes a unique vegetative structure and white cotton, may lend to this effect.
In comparison with the 645 nm spectral channel, lower black sky albedo values are found for both 470 nm and 555 nm channels for most crops. For example, maize field albedo values are around 0.1 for the 645 nm channel during the planting period while the values are around 0.075 and 0.05 for 555 nm and 470 nm spectral channels, respectively. Lower relative changes between planting and peak growing seasons are also found, with a 50% decrease in black-sky albedo for the blue 470 nm channel and a similarly muted 33% reduction in black-sky albedo in the green channel. This change is in line with expectations of albedo changes caused by photosynthesis activity as the majority of chloroplast absorption is found in the red and blue areas of the spectrum, with relatively less absorption of green by chlorophyll a and b and other absorptive pigments [18,19].
Plant matter has a high albedo at the 860 nm spectral channel and, therefore, as leaf area index (LAI) increases nearly opposite patterns from the visible channels are observed for the 860 nm channel (e.g., Figure 6). Lower black-sky albedo values on the order of 0.2 are found for maize, soybeans, and spring wheat during the planting season, while maximum mean black-sky albedo values around 0.4 are found during the peak growing season as identified in Figure 5. A similar, however more muted, increase is observed in black-sky albedo for cotton fields, with an increase from 0.25 to 0.35 from planting through peak growing season.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 18 the cotton plant, which includes a unique vegetative structure and white cotton, may lend to this effect.
In comparison with the 645 nm spectral channel, lower black sky albedo values are found for both 470 nm and 555 nm channels for most crops. For example, maize field albedo values are around 0.1 for the 645 nm channel during the planting period while the values are around 0.075 and 0.05 for 555 nm and 470 nm spectral channels, respectively. Lower relative changes between planting and peak growing seasons are also found, with a 50% decrease in black-sky albedo for the blue 470 nm channel and a similarly muted 33% reduction in black-sky albedo in the green channel. This change is in line with expectations of albedo changes caused by photosynthesis activity as the majority of chloroplast absorption is found in the red and blue areas of the spectrum, with relatively less absorption of green by chlorophyll a and b and other absorptive pigments [18,19].
Plant matter has a high albedo at the 860 nm spectral channel and, therefore, as leaf area index (LAI) increases nearly opposite patterns from the visible channels are observed for the 860 nm channel (e.g., Figure 6). Lower black-sky albedo values on the order of 0.2 are found for maize, soybeans, and spring wheat during the planting season, while maximum mean black-sky albedo values around 0.4 are found during the peak growing season as identified in Figure 5. A similar, however more muted, increase is observed in black-sky albedo for cotton fields, with an increase from 0.25 to 0.35 from planting through peak growing season. The black-sky albedo values at 1640 nm channel behave similarly to the visible channels but with a more drastic change in absolute magnitude between planting and peak growth periods. The changes in mean black-sky albedo values are from 0.3 to 0.2, 0.3 to 0.22, and 0.35 to 0.22 for maize, soybean, and spring wheat, respectively, from planting to peak growth periods. A 0.05 change in black sky albedo value from 0.35 to 0.3 is also found for cotton fields from the start of planting to peak growth. Similar patterns are also seen for the 2140 nm channel, although the absolute values are lower and the seasonal change is relatively more magnified. Surface albedo values of soil at these wavelengths is relatively high in comparison with healthy plants [20] which would explain the more muted response with respect to the lower wavelengths.
Distinct patterns, however, are observed for the 1240 nm spectral channels for most crops. In general, lower black-sky albedo values are found at the planting stage, with mean values of 0.25, 0.25, 0.3 and 0.33 for maize, soybean, spring wheat and cotton, respectively, that increases to approximately 0.35 for all four crop types around the peak of the growing season. Two local maximum black-sky albedo values at day 160 and day 220 are found for spring wheat, which represents a unique pattern relative to all other channels. Additionally, higher black-sky albedo The black-sky albedo values at 1640 nm channel behave similarly to the visible channels but with a more drastic change in absolute magnitude between planting and peak growth periods. The changes in mean black-sky albedo values are from 0.3 to 0.2, 0.3 to 0.22, and 0.35 to 0.22 for maize, soybean, and spring wheat, respectively, from planting to peak growth periods. A 0.05 change in black sky albedo value from 0.35 to 0.3 is also found for cotton fields from the start of planting to peak growth. Similar patterns are also seen for the 2140 nm channel, although the absolute values are lower and the seasonal change is relatively more magnified. Surface albedo values of soil at these wavelengths is relatively high in comparison with healthy plants [20] which would explain the more muted response with respect to the lower wavelengths.
Distinct patterns, however, are observed for the 1240 nm spectral channels for most crops. In general, lower black-sky albedo values are found at the planting stage, with mean values of 0.25, 0.25, 0.3 and 0.33 for maize, soybean, spring wheat and cotton, respectively, that increases to approximately 0.35 for all four crop types around the peak of the growing season. Two local maximum black-sky albedo values at day 160 and day 220 are found for spring wheat, which represents a unique pattern relative to all other channels. Additionally, higher black-sky albedo values are also found during and after harvesting seasons which is also in contradiction to patterns as observed from other channels; however, similar albedo values are found for planting and harvesting seasons.
The white-sky albedo values show similar behavior to the black-sky albedos, and similar figures as Figures 5 and 6 are included in the supplement documents as Figures S1 and S2. In summary, drastic changes in both black-and white-sky albedos are found, on the order of 30-75% for the major United States crops of maize, soybean, and spring wheat during the growth cycle for visible, NIR, and SWIR MODIS spectral channels. A similar, but lesser, change in magnitude is also found for cotton fields. Clearly, changes as significant as this in surface albedo during the plant growing period, as a function of plant type, may need to be taken into account when undertaking either remote sensing applications that require a knowledge of lower boundary conditions over cropland, or regional to global scale numerical weather/climate prediction and simulation models.
Of particular importance is the predictability of the annual cycle in cropland albedo as defined independently by each crop type. If there are significant annual inconsistencies in this cycle, on a per-crop basis, driven by external factors then the overall benefit of applying these albedo changes to future studies decreases substantially. To investigate the annual variations and their potential impacts, we separated each crop into its annual components and applied the earlier analysis over all wavelengths and albedo types independently. As observed in Figures 7 and 8 crops typically followed the same pattern of albedo changes, indicating the overall annual pattern as found from the study are recurring patterns.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 18 values are also found during and after harvesting seasons which is also in contradiction to patterns as observed from other channels; however, similar albedo values are found for planting and harvesting seasons. The white-sky albedo values show similar behavior to the black-sky albedos, and similar figures as Figures 5 and 6 are included in the supplement documents as figures S1 and S2. In summary, drastic changes in both black-and white-sky albedos are found, on the order of 30-75% for the major United States crops of maize, soybean, and spring wheat during the growth cycle for visible, NIR, and SWIR MODIS spectral channels. A similar, but lesser, change in magnitude is also found for cotton fields. Clearly, changes as significant as this in surface albedo during the plant growing period, as a function of plant type, may need to be taken into account when undertaking either remote sensing applications that require a knowledge of lower boundary conditions over cropland, or regional to global scale numerical weather/climate prediction and simulation models.
Of particular importance is the predictability of the annual cycle in cropland albedo as defined independently by each crop type. If there are significant annual inconsistencies in this cycle, on a percrop basis, driven by external factors then the overall benefit of applying these albedo changes to future studies decreases substantially. To investigate the annual variations and their potential impacts, we separated each crop into its annual components and applied the earlier analysis over all wavelengths and albedo types independently. As observed in Figures 7 and 8 crops typically followed the same pattern of albedo changes, indicating the overall annual pattern as found from the study are recurring patterns.

Variations in Broadband Shortwave (SW) Cropland Albedo Due to the Crop Growth Cycle
In addition to the individual channels, the MODIS BRDF product provides broadband albedo calculations in visible (0.3-0.7 nm), NIR (0.7-3.0 nm), and shortwave (SW; 0.3-3.0 nm) [9]. These have been compiled in a similar manner to the individual MODIS channels in this study with the mean and standard deviation calculated for each crop, HZ, and day of year. While the impact of crops on SW broadband albedo is less significant than the impacts on the individual channels, it is still nonnegligible. Figure 9a shows the mean broadband black-sky SW albedo for maize, and spring wheat for HZ4. Of the two crops, maize is found to have the larger seasonal swing in SW albedo, with an average SW black-sky albedo of 0.13 during the planting and early growing season with the SW black-sky albedo values changing to 0.17-19 during and after harvesting seasons. This still represents around a 30% change in broadband black-sky albedo during the crop-growing period, which could have a non-negligible impact on regional climates. A similar, but lesser, change in magnitude is found for spring wheat where the SW black-sky albedo values are observed to undergo a 10-20% change during the growing period. The relative change of SW black-sky albedo, however, is observed to be only ~10% for cotton fields. Very similar patterns are found for the white sky albedo, with an approximately 30% change in SW white-sky albedo for both maize and spring wheat as shown in Figure 9b.

Variations in Broadband Shortwave (SW) Cropland Albedo Due to the Crop Growth Cycle
In addition to the individual channels, the MODIS BRDF product provides broadband albedo calculations in visible (0.3-0.7 nm), NIR (0.7-3.0 nm), and shortwave (SW; 0.3-3.0 nm) [9]. These have been compiled in a similar manner to the individual MODIS channels in this study with the mean and standard deviation calculated for each crop, HZ, and day of year. While the impact of crops on SW broadband albedo is less significant than the impacts on the individual channels, it is still non-negligible. Figure 9a shows the mean broadband black-sky SW albedo for maize, and spring wheat for HZ4. Of the two crops, maize is found to have the larger seasonal swing in SW albedo, with an average SW black-sky albedo of 0.13 during the planting and early growing season with the SW black-sky albedo values changing to 0.17-19 during and after harvesting seasons. This still represents around a 30% change in broadband black-sky albedo during the crop-growing period, which could have a non-negligible impact on regional climates. A similar, but lesser, change in magnitude is found for spring wheat where the SW black-sky albedo values are observed to undergo a 10-20% change during the growing period. The relative change of SW black-sky albedo, however, is observed to be only~10% for cotton fields. Very similar patterns are found for the white sky albedo, with an approximately 30% change in SW white-sky albedo for both maize and spring wheat as shown in Figure 9b.
As suggested from Figure 9, variation in broadband albedo exists for different HZs and crops for a range of reasons including distinct planting date ranges, growing season length, and local crop varieties. For example, cropland planted in maize in HZ4 reaches the highest albedo around August 4th, whereas further south in HZ8 the peak in albedo can be found a full month earlier, around July 1st. Sensitivity to albedo timing and range can lead to significant changes in energy balance calculations, and accounting for these changes may improve long-term climate models [21,22]. As a result of this sensitivity, investigations into the long-term impacts of evolving cropping pattern effects on atmospheric changes should also consider effects of the overall albedo shift caused by the changes in the growth profile and reflectivity of each crop. Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 18 (a) (b) Figure 9. (a) Total shortwave, black-sky albedo by day of year for maize (red) and spring wheat (blue) in HZ4. (b) Total shortwave, white-sky albedo by day of year for maize (red) and spring wheat (blue) in HZ4.
As suggested from Figure 9, variation in broadband albedo exists for different HZs and crops for a range of reasons including distinct planting date ranges, growing season length, and local crop varieties. For example, cropland planted in maize in HZ4 reaches the highest albedo around August 4th, whereas further south in HZ8 the peak in albedo can be found a full month earlier, around July 1st. Sensitivity to albedo timing and range can lead to significant changes in energy balance calculations, and accounting for these changes may improve long-term climate models [21,22]. As a result of this sensitivity, investigations into the long-term impacts of evolving cropping pattern

Uncertainty Analysis
As mentioned in the dataset section, CDL data, which are constructed using satellite observations from satellites such as RESOURCESAT-1 and Landsat, contain uncertainties. Provided with the CDL dataset is the uncertainty of each classified crop type in the form of accuracy assessments. Based on these NASS accuracy assessments, the overall user accuracy of the CDL was calculated to be~80% on average for all land cover types each year; however, the major croplands utilized in this study generally fared better, with user accuracies averaging nearly 90% annually (as shown in popular US crops listed in Table 2). To minimize the impact of this uncertainty in crop type classification, only MODIS pixels that contained completely homogenous CDL croplands were selected for analysis. Areas that feature potential erroneous crop classification due to the difficulty of classification, observable as "noisy" spatial variations inside single field locations, such as field edges or variable wetlands, were also excluded. In addition, due to the extensive replication of both points and years, effects from erroneous classifications would be reduced in the final albedo calculations due to the sheer number of correctly classified croplands incorporated into the analysis. For the MODIS BRDF albedo product used in the study, evaluation efforts were undertaken to determine the validity of both the broadband and individual albedo channels under multiple conditions and testing regimes ( [24][25][26]), which found general agreement of other measurement techniques with the MODIS albedo. Still, it is possible that regions with the presence of thin clouds and aerosol layers may be misidentified as clear regions, as MODIS is insensitive to very optically thin cirrus clouds with optical depth less than 0.3 [27,28]. It is worth noting that the MODIS albedo product from the MODIS Multi-Angle Implementation of Atmospheric Correction (MAIAC) data [29] is also available. However, the MODIS MAIAC albedo product is only available every 8 days at a spatial resolution of 1 km. Yet, due to gridded infrastructure and property boundaries typical field sizes in the United States do not exceed 0.5 miles per side, effectively limiting the resolution to 800 m or finer for an albedo dataset to encompass MODIS pixels within single, homogeneous fields. Due to this spatial resolution limitation, as well as the reduced temporal resolution of 8 days between albedo calculations, the albedo product from the MAIAC data is not used.
Another area of concern is the level-3 nature of the MODIS BDRF product, which features a 16-day weighted method to estimate the albedo on the middle-most day of the sequence. This provides a continuous albedo curve for the entire reference period but can serve to smooth out both peaks and valleys of the albedo function, reducing accuracy at those transition points during the growing season. However, using level-3 gridded data offers the advantage of consistent size and location of the derived pixels, which allows the field locations to be constant throughout the study period, eliminating possible uncertainty arising from a daily changing pixel size and location as we would see with similar level-1, non-gridded, datasets.
Finally, while this study focused on crops in the United States, varieties and management practices vary greatly on a global scale. For example, planting date and irrigation practices of crops in regions that often double-crop (such as those seen in the southern Asia region [30]) can vary greatly from management practices found in the United States. Additionally, each individual crop can have several varieties, for example, as of February 2020 the USDA's Plant Variety Protection office lists 2134 varieties of field corn, 2752 varieties of soybean, and 1294 varieties of common wheat [31]. Due to the differences in both climate and management practices, these varieties can vary in growth speed, stage timing, and overall leaf area index, impacting the timing and scope of albedo changes throughout the season. As such, albedo changes derived from varieties found in the US may not be representative of albedo changes caused by varieties found outside the study area, even when factoring in planting dates, management, and weather conditions. As such studies of a specific region's crop growth patterns and timings should be completed to verify albedo patterns are similar to US results before incorporating the results from this database into studies outside this study's scope.

Evaluating the Feasibility of Using Changes in Normalized Difference Vegetation Index (NDVI) as a Proxy for Changes in Albedo over Cropland
In some of the previous studies, the changes in albedo over cropland during the growing seasons are approximated using NDVI values (e.g., Hsu, et al., [32]), particularly seen at the 650 nm and 860 nm wavelengths. The R 2 values between NDVI and black-sky albedo values at seven wavelengths (470 nm, 555 nm, 645 nm, 860 nm, 1240 nm, 1640 nm, and 2130 nm) are studied ( Figure 10) using four years' data as mentioned in Section 2. On a daily basis over MODIS land-cover defined cropland, a maximum R 2 value of 0.75 is found at 650 nm and the minimum R 2 value of 0.2 is found at 1240 nm. Also, as shown in Figure 10, this correlation on a daily scale increases over cropland from near zero correlation early in the growing season, reaching maximum correlation in the late-peak growth stage, and decreasing shortly thereafter back to near zero as crops are removed from the surface during harvest. eliminating possible uncertainty arising from a daily changing pixel size and location as we would see with similar level-1, non-gridded, datasets.
Finally, while this study focused on crops in the United States, varieties and management practices vary greatly on a global scale. For example, planting date and irrigation practices of crops in regions that often double-crop (such as those seen in the southern Asia region [30]) can vary greatly from management practices found in the United States. Additionally, each individual crop can have several varieties, for example, as of February 2020 the USDA's Plant Variety Protection office lists 2134 varieties of field corn, 2752 varieties of soybean, and 1294 varieties of common wheat [31]. Due to the differences in both climate and management practices, these varieties can vary in growth speed, stage timing, and overall leaf area index, impacting the timing and scope of albedo changes throughout the season. As such, albedo changes derived from varieties found in the US may not be representative of albedo changes caused by varieties found outside the study area, even when factoring in planting dates, management, and weather conditions. As such studies of a specific region's crop growth patterns and timings should be completed to verify albedo patterns are similar to US results before incorporating the results from this database into studies outside this study's scope.

Evaluating the Feasibility of Using Changes in Normalized Difference Vegetation Index (NDVI) as a Proxy for Changes in Albedo over Cropland
In some of the previous studies, the changes in albedo over cropland during the growing seasons are approximated using NDVI values (e.g., Hsu, et al., [32]), particularly seen at the 650 nm and 860 nm wavelengths. The R 2 values between NDVI and black-sky albedo values at seven wavelengths (470 nm, 555 nm, 645 nm, 860 nm, 1240 nm, 1640 nm, and 2130 nm) are studied ( Figure 10) using four years' data as mentioned in Section 2. On a daily basis over MODIS land-cover defined cropland, a maximum R 2 value of 0.75 is found at 650 nm and the minimum R 2 value of 0.2 is found at 1240 nm. Also, as shown in Figure 10, this correlation on a daily scale increases over cropland from near zero correlation early in the growing season, reaching maximum correlation in the late-peak growth stage, and decreasing shortly thereafter back to near zero as crops are removed from the surface during harvest. Improvements to this method can be found through limiting NDVI as a proxy to albedo to be considered only over homogeneous, non-barren crop fields. With this narrowing of scope gains in Improvements to this method can be found through limiting NDVI as a proxy to albedo to be considered only over homogeneous, non-barren crop fields. With this narrowing of scope gains in correlation are observed in the visible channels, with four-year R 2 values of 0.658, 0.504, and 0.661 found for the 0.47, 0.56, and 0.65 µm channels respectively; however, this improvement is not found in the NIR-SWIR range, with four-year R 2 values of 0.526, 0.030, 0.544, and 0.675 for the 0.86, 1.24, 1.64, and 2.13 µm channels respectively. These results can be further improved through implementation of the new database, which, in contrast to static or generic albedo databases, implements daily albedo changes for specific crops, incorporating regional management and growing patterns over homogenous croplands. Through the use of this database, most of the gains in the visible are maintained with four-year R 2 values 0.607, 0.536, and 0.621, and improves the correlation of the NIR-SWIR channels, with the exception of the 2.13 µm channel, with R 2 values of 0.650, 0.402, 0.620, and 0.686 for the 0.86, 1.24, 1.64, and 2.13 µm channels, respectively.
This exercise suggests that NDVI may not be an ideal proxy for capturing albedo changes over cropland during the growing seasons, particularly over mixed croplands and, as a result, further investigation into the impacts of surface crop type and land-use mixture on remote-sensing applications may be warranted.

Conclusions
Using four years (2015-2018) of CDL and MODIS BRDF data, the changes in black-and white-sky albedo at seven visible to shortwave IR wavelengths (470 nm, 555 nm, 645 nm, 860 nm, 1240 nm, 1640 nm, 2130 nm) as well as the broadband SW (300-5000 nm) during the planting, growing and harvesting periods were studied as functions of crop type and plant hardiness zones (HZ) over the continental United States. A total of 51 field crops were included in the study, with an emphasis on four commonly planted crops: maize, soybean, spring wheat and cotton. This study found: (1). Significant change in black-sky surface albedo at all seven channels are found between planting, peak growing, and harvest period for maize, soybean, and spring wheat. A total of a 70% decrease in 645 nm channel albedo is found for maize and soybean from early planting to the peak of the growing season, similar but less significant changes in the order of 60% and 30% are also found for spring wheat and cotton respectively. Additionally, similar but less significant changes are also found for the blue (470 nm) and green (555 nm) channels.
(2). For 1640 and 2140 nm spectral channels, although the overall black-sky albedo is much higher compared with the visible channels, lower albedo values are also observed at the peak growing period, with a~0.1 lower albedo found for the planting season for all four types of crop. The opposite pattern, however, is found for the 860 nm channels, with the peak surface albedo found during the peak growing season with a~0.2 albedo rise, compared to the planting season for maize, soybean, and spring wheat. A similar pattern, albeit weak in magnitude is found for cotton fields. Also note that the white-sky albedo behavior is nearly indistinguishable to the black sky albedo and therefore are not mentioned hereafter.
(3). An interesting pattern is found for the 1240 nm spectrum channel, with lower black-sky albedo values found for the planting seasons and higher black-sky albedo values found in both peak growing and harvesting seasons for maize, soybean, and spring wheat. The reason for this observed pattern is not fully known.
(4). In the broadband shortwave, both crop type and growth cycle were found to affect the overall black-sky albedo, varying by 30% relative (0.04 absolute) to the black-sky albedo at planting for maize, soybeans, and spring wheat, with cotton black-sky albedo changing only slightly, 10% relative (0.02 absolute). These significant variation in broadband albedo should be taken into consideration in climate applications over croplands in future studies. (5). We also found that NDVI may not be a good proxy for simulating cropland albedo changes during the growing cycle, as the R 2 values range from 0.2 to 0.75 between NDVI and black-sky surface albedo, for the seven spectral channels used in this study. This finding includes the key 555 nm channel which has an R 2 value of approximately 0.4 over mixed cropland.
Lastly, the goal of the study is to construct a cropland surface albedo database relating plant type, hardiness zone, and day of year to expected surface albedo using averages of nationwide, single-crop MODIS BRDF values over a four-year period (2015-2018) for a total of 51 field crops. This database provides the means to calculate appropriate albedo parameters for albedo-sensitive applications, such as long-term weather modeling, for the included cropland variations within the United States.
The database includes the mean, standard deviation, and pixel counts for matching crop, location, and time conditions for both white-and black-sky albedo allowing for a broad use of the data across multiple types of study. For example, linking the resulting albedo parameters to modeling of crop acreage changes driven by economic forces could allow estimation of currently unquantified surface albedo and regional climate effects of US agricultural and trade policies that affect farmers' choice of crops via the price mechanism (i.e., Starr et al. [33]). The database is attached as a supplement to this paper for potential remote-sensing, weather, and climate applications that require accurate quantification of lower boundary conditions. Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/18/2887/s1: Figure S1: Visible mean and standard deviation of WSA for four popular crops in selected HZs. Figure S2: IR mean and standard deviation of WSA for four popular crops in selected HZs. Albedo database.
Author Contributions: Authors J.Z. and J.S. designed the study. Author J.S. implemented the study. Authors J.S.R. and D.C.R. provided valuable comments during the study. All authors contributed to the writing of the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: This study is supported by the NSF project IIA-1355466.