Seasonal landcover estimate comparisons
We found that seasonal LULC estimates could have substantial differences between the growing and non-growing seasons. The 55 watersheds of our temperate study area had more tree cover during the growing season, typically by 5–10% of watershed area, according to the 2016 Dynamic World data (Fig. 2a). During the non-growing season, some areas classified as forested during the growing season were instead given built or shrubland LULC classes. Differences in seasonal LULC classifications in Dynamic World data were strongest in mixed-LULC watersheds (i.e., between 15% and 85% built LULC area), and weaker in very low built or very high built percentage watersheds (R2 = 0.50, df = 52, F = 25.64, p < 0.001; Fig. 2b). Built areas in particular can have substantial impacts on water quality due to the proliferation of impervious surfaces and increases in runoff and pollutants 36. Due to different classification approaches and resolutions, there is no exact comparison between 2016 Dynamic World built classification and the four NLCD 2016 developed classes. However, both datasets have somewhat similar estimates of developed land. There is a relative mean absolute difference (RMAD) in built land of 9.4% of watershed area between NLCD 2016 and Dynamic World 2016 growing season data, and 6.0% between NLCD 2016 and Dynamic World 2016 non-growing season data (Fig. S1 and Table S1 in the Supplementary Material).
Our study watershed for Case #2 (the Bush Creek watershed of Maryland, USA) showed an 8% increase in built LULC, and a 16% decrease in tree cover, in non-growing season relative to growing season Dynamic World data from 2016. A similar seasonal cycle in built LULC classification was observed over the entire time period of available Dynamic World data (Fig. 2c,d,e). This confirms that the differences between seasons were not limited to a single year of data. Meanwhile, our study watershed for Case #3 (the Difficult Run watershed of Virginia, USA) showed a 12% decrease in tree cover and a 10% increase in built areas in the non-growing season compared to the growing season Dynamic World data from 2016. Also, some watersheds had LULC differences in crop and grassland cover between seasons, usually less than 10%, suggesting there could be changes in agricultural watersheds in addition to the more urbanized and forested watersheds of this study (Fig. 2a).
LULC changes between seasons were often concentrated along forested edges of mixed-LULC areas (Fig. S2 in the Supplementary Material). In these deciduous areas, such as the edges of mixed residential/forested zones in Fig. S2, leaf cover decreases during the non-growing season, which could be exposing other types of LULC underneath, or making forest more difficult to distinguish from surrounding built area. In prior studies, characteristics such as atmospheric conditions (e.g., aerosol scattering, water vapor, absorption of light) and reflectance (e.g., bidirectional reflectance and zenith angle) have also been found to cause non-random errors in top-of-atmosphere readings. These errors could cause differences in spectral data compared to corrected surface reflectance products for classifying landcover such as vegetation 29,30,37. Dynamic World used a calibrated surface reflectance product to train the classifier (Sentinel-2 Level-2A; L2A) and a top-of-atmosphere product (Sentinel-2 Level 1C; L1C) to generate the dataset 4. L2A was chosen to train the classifier because it accounts for characteristics like solar distance, zenith, and atmospheric conditions, but could not be used to generate the LULC dataset because it does not cover the entire time period of Dynamic World data 4. Thus, any seasonal cycles in solar distance, zenith, and atmospheric conditions would not be corrected for in the L1C product used to classify Dynamic World data 4.
Case #1: Water quality regression models
We found that regression-based water quality models for specific conductance (SC) fit to different seasons of LULC data could show similar relationships. We focused on SC of the water as a dependent variable since it is one of the most accessible and widely used indicators of water quality. Median stream water SC was positively correlated with 2016 Dynamic World built LULC during all seasons (Fig. 3a), as would be expected due to the influence of point- and non-point sources of salt pollution from developed areas 36.
Regression models comparing growing and non-growing season Dynamic World 2016 LULC data with these median SC values over the entire time period did not show large differences between seasons, as shown in Table S2 of the Supplementary Material and Fig. 3a. For example, the model for growing season built LULC vs. median SC had a slope of 6.39, while the same model for non-growing season LULC had a slope of 6.06, and the AIC’s (metrics of performance where higher AIC’s are better) for both models were within 2 AIC units, which suggests similar performance. For perspective, a model created with developed classes from NLCD 2016 had a slope of 6.23 (Table S2).
This demonstrates regression models fit with sub-annual Dynamic World LULC data could comparably estimate general water quality tendencies in watersheds. This would be expected because processes such as salinization of surface waters often operate at longer temporal scales 38,39. For heavily urbanized areas with high built percentage, however, the relationship between built LULC and specific conductance could be less clear, as Dynamic World pixels do not distinguish between types of built developments, for instance between single-family housing and commercial or industrial land. Our approach is limited in that the National Park Service water quality dataset is based on discrete sampling of relatively low frequency (monthly to every three months). Benefits or limitations of the higher spatial and temporal resolution of the Dynamic World data might become more apparent if exploring sensor-derived high-frequency water quality datasets, which are growing rapidly across the US and worldwide, particularly for SC.
Case #2: Hydrologic/water quality models for ungauged basin
We found that there could be large differences in the outputs of uncalibrated hydrologic models that demonstrate objectively how seasonal differences in LULC estimates could be problematic for interpretation of results. The uncalibrated (meaning not optimized to a specific hydrological system) SWAT model for the ungauged Bush Creek watershed simulated 19% greater surface runoff (Fig. 3b), 42% greater soluble phosphorus yield (Fig. 3c), and 51% greater nitrate runoff (Fig. 3d) over the entire 2005–2018 period when non-growing season Dynamic World 2016 data was used as the LULC input, rather than growing season LULC. These differences are relatively greater than the current pollutant load reduction targets for Chesapeake Bay of 17% total nitrogen load and 12% phosphorus load 40, highlighting the importance of considering the dynamics of seasonal LULC in order to effectively design water quality improvement efforts, particularly when decision making is involved.
The observed differences are due to the 8% increase in built areas in non-growing season data, which have more impervious surfaces, a higher runoff curve number, and generate proportionally more runoff than the forested areas which were classified during the growing season. This demonstrates that in watersheds with mixed-LULC areas, seasonal differences in LULC classification could result in substantial changes to model outputs for ungauged or uncalibrated basins. This could be particularly problematic when using coarse-scale, computationally efficient SWAT models that assign subbasin conditions based on the dominant hydrologic response unit, as a change in dominant LULC type in a watershed would result in different subbasin conditions in the model greater than the proportional change in LULC. For perspective, these outputs were lower than those simulated by the uncalibrated SWAT model with NLCD 2016 LULC input (Fig. 3b,c,d). Nitrate and soluble phosphorus yields are lower in the uncalibrated Bush Creek watershed models with Dynamic World 2016 inputs than in the uncalibrated model with NLCD 2016 input because the NLCD 2016 classification included more developed and agricultural lands, which increase SWAT nutrient yields (Table S1 in the Supplementary Material). NLCD 2016 also provides a greater breakdown of LULC classes. For instance, the NLCD 2016 open space class (18% of watershed area) was treated as developed area in the SWAT model and thus assigned elevated nutrient concentrations. This resulted in 10 of the 31 subbasins having urban as the dominant LULC in the model with NLCD 2016, compared with 6 urban-dominant subbasins in the model with non-growing season Dynamic World 2016, and only 2 urban-dominant subbasins in the model with growing season Dynamic World 2016. Also, the three models with different LULC inputs simulated differences in average annual water yields (438 mm with growing season Dynamic World 2016, 466 mm with non-growing season Dynamic World 2016, and 455 mm with NLCD 2016), which would impact nutrient yields.
This supports the previously stated need for improved seasonal process representation for water quality and hydrologic models such as SWAT, particularly for those used in ungauged catchments that have no stream measurement stations to compare to 41. This information could be useful for informing urban planners about the impacts of LULC change on flash flood disasters in ungauged watersheds using hydrologic models such as SWAT 42. This also could be applicable for other common environmental models that are used to simulate water quality/quantity in ungauged or uncalibrated basins and often use growing-season LULC datasets, such as the Generalized Watershed Loading Functions (GWLF) model and Model My Watershed app 43–45 or models of the Continental Hydrologic Intercomparison Project, which uses the PARFLOW-CONUS model and National Oceanic and Atmospheric Administration National Water Model configuration of the WRF-HYDRO model 46–48. A brief summary of the performance of uncalibrated models for simulating observed streamflow data can be found alongside the calibrated model results in Case #3.
Case #3: Calibrated hydrologic models
We found that optimization of hydrologic parameters for models using LULC data from different seasons could result in comparable performance between the models for simulating an observed dataset, however the parameter values would differ depending on the LULC season calibrated for. The individually calibrated SWAT models (meaning they had been optimized to best represent measured streamflow at the watershed outlet independently of each other) using growing season vs. non-growing season Dynamic World 2016 LULC input for the Difficult Run watershed had comparable performance for simulating streamflow, despite the differences in LULC inputs (10% increase in built areas and 12% decrease in tree cover for the non-growing season LULC input). For instance, Nash Sutcliffe Efficiency (NSE) performance metrics at the daily time step were between 0.52 and 0.54 for each model with Dynamic World LULC over the calibration and validation time periods, while Kling-Gupta Efficiency (KGE) was between 0.61 and 0.75, and refined Index of Agreement (dr; which by not squaring errors provides a better measure of low flow performance) only ranged between 0.68 and 0.70 (Table 1). Higher values for these metrics suggest better performance. These values are similar to those of previous SWAT evaluations in urban watersheds that occurred at monthly time steps 49,50, and at the daily time step our models had to undergo much higher temporal scrutiny. For perspective, the SWAT model calibrated with NLCD 2016 LULC input had an NSE of 0.50 for the calibration period and 0.48 over the validation period, which was lower performance than the models with Dynamic World 2016 growing or non-growing season LULC data in this case (Table 1), and could reflect difficulties simulating hydrology in urban areas with complex stormwater pathways, as the Difficult Run watershed was 58% developed area in the NLCD 2016 data. Visually, the calibrated models with Dynamic World 2016 and NLCD 2016 LULC performed similarly to each other, although the non-growing season Dynamic World data appears to have done a better job of simulating low streamflow events in this case, which could contribute to the slightly higher dr values (Fig. 3e,f). The reason for the comparable performance was that the model parameters (i.e., coefficients in equations) adjusted to account for the change in LULC as it relates to the streamflow calibration, as has been previously documented with hydrologic models such as SWAT 33.
Table 1
Comparison of performance for calibrated SWAT model developed with Dynamic World 2016 growing season landcover input, calibrated SWAT model developed with Dynamic World 2016 non-growing season landcover input, and calibrated SWAT model developed with NLCD 2016, over the calibration (2005–2011) and validation (2011–2018) periods, at the daily time step for the Difficult Run watershed. Also includes performance of uncalibrated models over the validation period. Performance indices are R2, Nash Sutcliffe Efficiency (NSE), Kling-Gupta Efficiency (KGE), and refined Index of Agreement (dr).
SWAT landuse input | Period | R2 | NSE | KGE | dr |
Growing season | Uncalibrated | 0.39 | -1.04 | 0.05 | 0.52 |
Non-growing season | Uncalibrated | 0.40 | -1.08 | 0.03 | 0.50 |
NLCD 2016 | Uncalibrated | 0.36 | -0.93 | 0.15 | 0.52 |
Growing season | Calibration | 0.54 | 0.53 | 0.61 | 0.69 |
Non-growing season | Calibration | 0.54 | 0.54 | 0.65 | 0.70 |
NLCD 2016 | Calibration | 0.51 | 0.50 | 0.55 | 0.69 |
Growing season | Validation | 0.56 | 0.53 | 0.73 | 0.68 |
Non-growing season | Validation | 0.57 | 0.52 | 0.75 | 0.68 |
NLCD 2016 | Validation | 0.53 | 0.48 | 0.67 | 0.67 |
This emphasizes the importance of having a mechanistic understanding of how optimized model parameters are dependent on LULC scenarios and that LULC changes can impact simulation of hydrological processes 33. Based on PAWN sensitivity analysis results, the most sensitive parameters for our Difficult Run study were channel hydraulic conductivity (CH_KII), bank flow recession coefficient (ALPHA_BNK), and runoff curve number (CN_F) (Figure S3 in the Supplementary Material). Among these and other sensitive parameters, there were differences in optimized values depending upon the SWAT LULC input. For example, ALPHA_BNK was optimized to 0.14 for growing season Dynamic World 2016, 0.27 for non-growing season Dynamic World 2016, and 1.00 for NLCD 2016 data (Table S3 of the Supplementary Material). Also, CN_F adjustment optimized to -0.17, -0.20, and − 0.12, respectively, for these three data inputs, showing that the optimization adjusted runoff processes to compensate for the different proportions of LULC types in the watershed. Previous work has also found that SWAT model optimized parameters can be sensitive to different LULC scenarios, in particular the conversion of forests to other types of LULC 33. There was a difference in forests of 12% of watershed area between growing and non-growing season Dynamic World 2016 data for the Difficult Run watershed in our study (Table S1 in the Supplementary Material), which is as large a difference as real changes in forests that have been found to cause these sensitivities in model parameters 33. Thus, the seasonal difference in forest type LULC estimates is one possible explanation for our differences in optimized parameter values.
The calibrated Difficult Run watershed models had better performance than the uncalibrated models, which for instance had an NSE of -1.04 for growing season Dynamic World 2016 LULC input and − 1.08 for non-growing season LULC input during the validation period (a separate time period of data to evaluate models against more objectively) (Table 1). Prior to calibration, the SWAT model with non-growing season LULC input had particularly poor performance, likely because the default parameter values did not account for the inherent increase in built areas or decrease in forests observed in non-growing season LULC data. Thus, we did not use the uncalibrated models beyond objectively demonstrating the SWAT model’s sensitivity to LULC.
Also, although the calibrated models with different Dynamic World LULC season inputs achieved comparable performance for simulating streamflow over the validation time period, it is important to consider that the differences in parameter values creates the potential for the models to respond differently to future changes in LULC or climate change due to variations in unmeasured water balance outputs 31. An approach such as time-varying parameters may help accommodate these changes, particularly over long time periods 33. Examples of other calibrated environmental models these findings could be applicable to are the 13 models of the Great Lakes Runoff Intercomparison Study 51, models of water temperature based on hydrologic fluxes 52, and studies using other calibrated geospatial models such as SPARROW 53,54 and ParFlow 55,56.