The New England region of the northeastern United States has a land use history characterized by forest clearing for agriculture and other uses during European colonization and subsequent reforestation following widespread farm abandonment. Despite these broad changes, the potential influence on local and regional climate has received relatively little attention. This study investigated wintertime (December through March) climate impacts of reforestation in New England using a high-resolution (4 km) multiphysics ensemble of the Weather Research and Forecasting Model. In general, the conversion from mid-1800s cropland/grassland to forest led to warming, but results were sensitive to physics parameterizations. The 2-m maximum temperature (T2max) was most sensitive to choice of land surface model, 2-m minimum temperature (T2min) was sensitive to radiation scheme, and all ensemble members simulated precipitation poorly. Reforestation experiments suggest that conversion of mid-1800s cropland/grassland to present-day forest warmed T2max +0.5 to +3 K, with weaker warming during a warm, dry winter compared to a cold, snowy winter. Warmer T2max over forests was primarily the result of increased absorbed shortwave radiation and increased sensible heat flux compared to cropland/grassland. At night, T2min warmed +0.2 to +1.5 K where deciduous broadleaf forest replaced cropland/grassland, a result of decreased ground heat flux. By contrast, T2min of evergreen needleleaf forest cooled –0.5 to –2.1 K, primarily owing to increased ground heat flux and decreased sensible heat flux.
Changes in forest cover affect climate by altering biogeophysical surface properties that include both radiative (e.g., albedo) and nonradiative (e.g., surface roughness and evapotranspiration) forcings (Lee et al. 2011; Davin and de Noblet-Ducoudré 2010; Bonan 2008). The net contribution of these forcings is difficult to quantify in the midlatitudes for several reasons, including seasonal variations in climate response (Bonan 2008) and nonlinear effects that may result from interactions among forcings (Davin and de Noblet-Ducoudré 2010). Wintertime responses are of particular interest in the midlatitudes due to the high albedo of snow (0.7–0.8) over deforested lands in comparison to the low albedo of snow-covered forests (0.2–0.3; Jin et al. 2002; Betts and Ball 1997; Robinson and Kukla 1984).
Surface albedo, the ratio of reflected to incoming solar radiation, generally increases as a result of deforestation in northern midlatitude winter due to the removal of low albedo forest canopies that would otherwise mask highly reflective snow cover (Betts and Ball 1997; Betts et al. 2007; Kvalevåg et al. 2010; Davin and de Noblet-Ducoudré 2010; Robinson and Kukla 1985). Numerous global and continental modeling studies have identified albedo as the main driver of high and midlatitude cooling in response to deforestation (Betts 2001; Betts et al. 2007; Feddema et al. 2005; Pitman et al. 2009; Davin and de Noblet-Ducoudré 2010; Kvalevåg et al. 2010). Betts (2001) reported −1 to −2 K of cooling in winter and spring in northern midlatitudes due to the increase in surface albedo over deforested lands covered with snow. Similarly, Feddema et al. (2005) reported regional cooling up to −2 K due to increases in surface albedo by replacing midlatitude forests with cropland. Most recently, the “Land-Use and Climate, Identification of Robust Impacts” (LUCID) project compared seven global climate models and reported an interquartile range of −0.2 to −0.9 K of cooling in North America winter (December–February) due to an increase in crop and pasture between 1870 and 1992 (de Noblet-Ducoudré et al. 2012).
The nonradiative impacts of deforestation in winter are less well known. The reduction in surface roughness from deforestation generally decreases turbulence in the boundary layer, although complex nonlinear interactions within the boundary layer make it challenging to model this phenomenon (Fernando and Weil 2010; Delage et al. 2002). The simulation of a stable boundary layer (e.g., inversions) over cold surfaces such as snow and sea ice remains an active area of research (Delage 1997; Derbyshire 1999; Sterk et al. 2013), yet this nonradiative process could contribute more to the observed cooling associated with midlatitude deforestation than surface albedo (Davin and de Noblet-Ducoudré 2010; Lee et al. 2011; Boisier et al. 2012; Luyssaert et al. 2014; Zhao and Jackson 2014).
Diurnal temperature responses were evaluated in a study comparing above-canopy tower measurements to nearby climate stations located in open lands (Lee et al. 2011). For sites located between 28° and 45°N, January mean temperatures were +0.52 K warmer over forests, with stronger warming at night compared to day. During the day, warming due to radiative forcing of low albedo over forest is compensated by cooling related to energy redistribution from higher surface roughness over forest, where enhanced mixing during the day causes surface temperature of forests to more efficiently dissipate heat compared to open lands that suppress mixing. At night, stable conditions over open lands with low surface roughness trap cold air at the surface; the higher surface roughness of forests increases turbulence and draws heat toward the forest canopy from aloft (Lee et al. 2011).
Few studies have focused on the climate effects of reforestation in the northeastern United States. While much of the midlatitudes have remained deforested to the present day (Pitman et al. 2009), the forests in the northeastern United States have followed a recovery trajectory. When early European settlers first arrived in the United States in the 1600s, 90% of the New England region in the northeastern United States (Fig. 1a) was covered with forest (Fig. 1b; Foster et al. 2008). By 1850, forest-covered area had decreased to just over 50%, with forest having been removed for pasture, agriculture, lumber and timber, and fuel (Foster et al. 2008). In several states (e.g., Connecticut, Rhode Island, and Massachusetts) forest cover dropped to below 35%. Farmland abandonment and alternative fuel sources (e.g., coal) led to extensive reforestation in the New England region through the mid-1900s (Baldwin 1942; Foster et al. 2008). Recently, forest cover in New England has again begun to decline, primarily due to urban and suburban development and mechanical disturbance (e.g., clear cutting; Barnes and Roy 2008). Currently, New England forest cover is around 75% (Fig. 1b).
Robinson and Kukla (1984) estimated a doubling of surface albedo as a result of mid-1800s deforestation in the central and eastern United States, but they did not assess the associated impacts on surface temperature. Klingaman et al. (2008) used the fifth-generation Pennsylvania State University–National Center for Atmospheric Research Mesoscale Model (MM5) in central and northern Pennsylvania, replacing presettlement evergreen needleleaf forest with barren and sparse land cover to evaluate climate responses to deforestation in five 1-month February (2000–04) simulations. Results suggested the responses in winter temperature were small and variable, possibly due to the relatively shallow snowpack during the five February simulations. The mean difference in albedo between the Klingaman et al. (2008) forested and deforested scenarios was 0.15, whereas observational studies indicate that the difference between forest and open land snow-covered albedo is closer to 0.45 (e.g., Jin et al. 2002; Betts and Ball 1997; Robinson and Kukla 1984). Additionally, the simulations may not have captured the nighttime cooling over open land relative to forest that is associated with changes in surface roughness proposed by Lee et al. (2011). Furthermore, only one land surface model (LSM), the Noah LSM, was used in the Klingaman et al. (2008) study. The magnitude and detection of climate responses to land cover change can be very sensitive to the underlying surface datasets and parameterization in land surface models (Oleson et al. 2004).
In this study, we simulated cold season (November through April) climate using the Advanced Research version of the Weather Research and Forecasting (WRF-ARW; Skamarock et al. 2008) Model and compared a present-day land cover scenario to a mid-1800s deforested land cover scenario. Simulations included a 12-member physics ensemble that used three land surface models, two radiation schemes, and two microphysics schemes. We evaluated the robustness of responses in winter (December through March) minimum temperature, maximum temperature, and surface energy fluxes due to forest cover change for a historically cold, snowy winter (December 2008–March 2009) and a warm, dry winter (December 2011–March 2012).
2. Methods and datasets
a. Weather Research and Forecasting Model
We used the WRF-ARW (Skamarock et al. 2008) to simulate winter climate over the eastern United States (Fig. 2). The modeling domain included triple one-way nests over the following three domains: d01: 36-km grid spacing over the eastern United States, d02: 12-km grid spacing over the northeastern United States, and d03: 4-km grid spacing over the New England region.
Lateral boundary conditions came from the ERA-Interim dataset (Dee et al. 2011) and provided the large-scale atmospheric forcing to drive the WRF Model. The simulations included a 35-day spinup (27 October–30 November) that preceded the analysis and was necessary to initialize soil temperatures for a snow-free surface. We simulated climate responses to land cover change for two climatic extremes within the past decade, a “cold, snowy” winter (December 2008–March 2009) and a “warm, dry” winter (December 2011–March 2012) (Table 1; Northeast Regional Climate Center; http://www.nrcc.cornell.edu). The purpose of using two extreme seasons was to identify whether the simulation of winter extremes (e.g., high vs low snow cover) is sensitive to physics parameterizations. Greenhouse gas concentrations were held at present-day values in order to isolate the effects of historical land cover change on surface energy fluxes.
b. Land cover scenarios
Winter climate was simulated using two land cover scenarios: 1) present-day United States Geological Survey (USGS) 24-class land cover (Fig. 2b) and 2) historical 1850 deforested land cover (Fig. 2c). The historical mid-1800s deforested scenario is a modified version of the History Database of the Global Environment version 3.1 (HYDE3.1) historical land use data for the year 1850 (Klein Goldewijk et al. 2011). The HYDE3.1 dataset provides fractional cropland and grassland area at 5-min resolution. It uses improved subnational data to map land cover based on urban and rural population and implements allocation algorithms with time-dependent weighting maps for cropland and grassland. HYDE3.1 historical 1850 land cover data, however, underestimates mid-1800s deforestation in northern New England compared to more robust estimates of forest cover for the region that use 1850 census data of improved (e.g., cropland and pasture) and unimproved (e.g., woodlots) acreage (Foster et al. 2008; Baldwin 1942; Harper 1918). We describe below how we modified HYDE3.1 to increase the 1850 percentage deforested area to reflect the forest cover area used by Foster and coauthors (Foster and Aber 2004; Foster et al. 2008, 2010; Fig. 1b).
The deforestation modification process first identifies present-day USGS 24-class land cover grid cells as candidates for mid-1800s deforestation using HYDE3.1. Next, starting from the lowest elevation, non-water-body grid cells were converted from present-day land cover to cropland/grassland mosaic increasing in elevation until the percentage forest cover within a given New England state boundary matched the percentage forested area presented in Fig. 1b (Foster et al. 2008). In the event of a tie in elevation, present-day urban and built-up and forested land cover classes were preferentially converted to cropland/grassland over forested areas.
c. Multiphysics ensemble configurations
We use a multiphysics ensemble to explore model uncertainty in the climate responses to the mid-1800s and present-day land cover scenarios. Twelve physics configurations (Table 2) were run for each of the two climate extremes and two land cover scenarios, for a total of 48 simulations. Physics options included two longwave/shortwave (LWSW) schemes, two microphysics (MP) options, and three land surface model schemes, described in more detail below.
The two microphysics schemes included the WRF single-moment six-class scheme (MP6; Hong and Lim 2006) and the Thompson graupel scheme (MP8; Thompson et al. 2008). The two schemes are similar in the number and types of hydrometeors they simulate (rain, hail, snow, water vapor, cloud ice, and graupel) but differ in assumptions of snow size distributions and snow shape (spherical vs nonspherical).
The first LWSW combination evaluated was the Rapid Radiative Transfer Model (RRTM; Mlawer et al. 1997) longwave scheme, which uses look-up tables for absorption at multiple bands for CO2 (379 ppm), N2O (319 ppb), and CH4 (1774 ppb), with the Goddard shortwave scheme, a two-stream multiband scheme with ozone from climatology and cloud effects (Chou and Suarez 1994). This LWSW combination is hereafter referred to as RRTMG. The second LWSW scheme (here after CAM3.1) included the Community Atmosphere Model (CAM3.1) longwave scheme, used with the CAM3.1 shortwave scheme. The CAM3.1 longwave scheme accounts for aerosols and trace gases, annual CO2 changes, and constant N2O (311 ppb) and CH4 (1714 ppb) concentrations. The CAM3.1 shortwave scheme uses δ-Eddington approximation to simulate effects of multiple scattering (Collins et al. 2006). The LWSW combinations chosen here are most suitable for regional climate simulations (Mooney et al. 2013).
The three land surface models included 1) the Community Land Model Version 4.0–Satellite Phenology (CLM4.0-SP; Oleson et al. 2010), 2) the Noah Multiple Parameterization (Noah-MP; Niu et al. 2011) with Biosphere Atmosphere Transfer Scheme albedo (NoahMP-BATS), and 3) the Noah-MP with Canadian Land Surface Scheme albedo (NoahMP-CLASS).
CLM4.0 is called as a subroutine in WRF3.5.1 (Zhao et al. 2014; Lu et al. 2015). The 24-class USGS land cover types are converted into one of five subgrid land cover types—lake, wetland, urban, glacier, or vegetated. Vegetated subgrid land cover types are subsequently assigned to up to four of 16 available plant functional types (PFTs). Monthly leaf area index (LAI) and stem area index (SAI) are prescribed for each PFT. When snow is present, LAI and SAI are adjusted for the vertical fraction of vegetation buried by snow, and 0.2 m is the snow depth at which short vegetation is considered fully buried by snow (Oleson et al. 2010). Snow albedo is simulated with the Snow, Ice, and Aerosol Radiative Model (SNICAR; Flanner and Zender 2006; Flanner et al. 2007), which uses a two-stream radiative transfer solution (Toon et al. 1989) and includes effects of snow grain size, solar zenith angle, and interannually varying snowpack impurities such as black carbon, soot, and aerosols.
The Noah-MP builds on the original Noah LSM by including multiple options to parameterize up to 12 surface related processes, summarized in Table 3. The two albedo schemes tested include the 1) Biosphere–Atmosphere Transfer Scheme (BATS; Yang et al. 1997) albedo parameterization and 2) the Canadian Land Surface Scheme (CLASS; Verseghy 1991) albedo parameterization. The BATS scheme is a sophisticated scheme that includes calculation of snow albedo for direct and diffuse radiation over visible and near-infrared wave bands, and accounts for the effects of solar zenith angle, grain size growth, and snowpack impurities such as dirt or soot on snow on snow age. The CLASS scheme is a simpler computation that accounts for the decrease in snow albedo as a snowpack ages. Of the two Noah-MP albedo options tested, the BATS scheme is most similar to CLM4.0.
d. Model validation data sources
1) PRISM monthly temperature and precipitation
Model skill was evaluated for simulations using the present-day land cover scenario. The Parameter-Elevation Relationship on Independent Slopes Model (PRISM) daily maximum temperature, minimum temperature, and precipitation data were chosen as a benchmark for model validation (Daly et al. 2008; Di Luzio et al. 2008). PRISM uses the National Elevation Database 3-arcsecond digital elevation model (DEM) to generate climate–elevation regression for 13 000 precipitation and 11 000 temperature stations that are assigned weights based on factors including elevation, proximity to the coast, topographic slope and aspect, and orographic effectiveness of the terrain. PRISM data were regridded from their native 30-arcsecond resolution (~800 m) to the 4-km modeling domain encompassing the New England states (Fig. 2b) using bilinear interpolation for temperature and first-order conservative interpolation for precipitation.
PRISM is commonly used for model evaluation in the United States (e.g., Chen et al. 2014; Lu and Kueppers 2012; Rasmussen et al. 2011; Subin et al. 2011). PRISM is based on surface observations from climate stations that typically measure temperature “at a height of between 1.2 m to 2 m above ground level… and not shielded by, or close to, trees …” (WMO 2008, p. 1.2-3). PRISM temperature is therefore primarily representative of deforested landscapes and does not necessarily represent temperature over forested lands. The WRF Model provides 2-m minimum (T2min) and maximum (T2max) temperature as output fields that are compared to PRISM.
2) NOHRSC/SNODAS snow depth data
The cold, snowy and warm, dry modeled snow depth (SNOWH) data from the present-day land cover WRF ensemble members were evaluated against monthly 1-km gridded snow depth data from the NOAA National Weather Service’s National Operational Hydrologic Remote Sensing Center (NOHRSC) Snow Data Assimilation System (SNODAS) dataset (National Operational Hydrologic Remote Sensing Center 2004; Carroll et al. 2001). SNODAS is a modeling and data assimilation system developed by NOHRSC to provide the best possible estimates of snow parameters. SNODAS integrates satellite-, airborne-, and ground-based snow data with model estimates of snow cover (Carroll et al. 2001). The 1-km gridded SNODAS snow depth data were regridded using bilinear interpolation to 4-km for comparison with WRF Model output.
3) Surface albedo validation data
We compared WRF modeled albedo to the 500-m Moderate Resolution Imaging Spectroradiometer (MODIS) BRDF/Albedo product version 005 (MCD43A3; Schaaf et al. 2002). MODIS snow-covered and snow-free shortwave broadband (0.3–5.0 μm) albedo statistics (mean, standard deviation) were generated for five MODIS Land Cover version 005 (MCD12Q1; Friedl et al. 2010) International Geosphere Biosphere Programme (IGBP) classes: 1) cropland, 2) grassland, 3) mixed forest, 4) deciduous broadleaf forest, and 5) evergreen needleleaf forest for comparison with WRF modeled albedo. MODIS cropland and grassland albedo were averaged for comparison with the WRF cropland/grassland mosaic land cover class.
Daily ground-based observations of snow-covered and snow-free albedo were collected by the Community Collaborative, Rain, Albedo, Hail, and Snow (CoCoRAHS) volunteer network in New Hampshire (Burakowski et al. 2013; data available at www.cocorahs-albedo.org). CoCoRAHS daily snow depth and albedo data collected between December 2011 and March 2014 were used to evaluate the relationship between albedo and snow depth in the three land surface models tested.
a. Model validation
The present-day land cover ensembles were evaluated to assess how well WRF3.5.1 performs relative to PRISM data. Comparisons between WRF and PRISM seasonal T2max and T2min are summarized in a Taylor diagram (Fig. 3; Taylor 2001). WRF simulations using the CLM4.0 and RRTMG radiation schemes produced the highest correlations and lowest normalized standard deviations in T2max for both winters. The dominant factor contributing to biases in modeled T2max was the choice of land surface model and, to a lesser extent, radiation scheme (Figs. 3a,b). The microphysics scheme did not influence T2max biases. Considerable warm biases in T2max were observed in the NoahMP-BATS and NoahMP-CLASS simulations for both the cold, snowy winter (Fig. S1a) and warm, dry winter (Fig. S1b). Warm biases in NoahMP were greatest (+6 to +8 K) in the northeastern part of the domain, where mixed forest and evergreen needleleaf forest are the primary land cover types (Fig. S1). All ensemble members had a cool bias relative to PRISM for both the cold, snowy and warm, dry winters over present-day agricultural land in northeastern Maine (Fig. S1). Cold biases in T2max in the CLM4.0 simulations were stronger in simulations that used the CAM3.1 radiation scheme, reaching −4 K in the southwestern part of the domain and strongest in the cold, snowy simulation (Fig. S1a). Biases in T2max were relatively weak (±2 K) in the CLM4.0 simulations run with the RRTMG radiation scheme.
Daily T2min biases were influenced to a greater extent by the LWSW scheme than by the land surface model (Figs. 3c,d). Each of the 12 physics configurations produced reasonable simulations of winter T2min, although warm biases were greater in the warm, dry simulations for ensemble members using RRTMG (Fig. S2). The CAM3.1 scheme produced cool biases up to −5 K in the southern part of the domain in the cold, snowy simulation (Fig. S2a) and in the northern part of the domain in the warm, dry simulation (Fig. S2b).
None of the multiphysics ensemble members adequately simulated total winter precipitation (Fig. 4). Seasonal (December–March) precipitation biases were slightly higher using the MP6 microphysics scheme compared to the MP8 scheme (Fig. S3). The strongest positive biases in total precipitation tended to occur over regions of steep terrain. Specifically, the Green Mountains in Vermont, the White Mountains in New Hampshire, and Mt. Katahdin in north central Maine consistently overpredicted precipitation (>50%) in nearly all simulations (Fig. S3). The combination of MP8 microphysics and the CAM3.1 schemes produced lower precipitation biases in steep terrain in both the cold, snowy and warm, dry simulations.
3) Snow depth
Modeled snow depth was dependent predominantly on the choice of LWSW scheme. The CAM3.1 scheme produced greater positive snow depth biases than the RRTMG scheme in the southern part of the domain during the cold, snowy simulation, up to 0.4 m (Fig. S4). Negative snow depth biases relative to SNODAS tended to occur along the northern border of Maine in all ensemble members.
Snow-covered and snow-free albedo averages were calculated for the three different land surface models. Snow-covered albedo over cropland/grassland mosaic grid cells was 0.73 ± 0.15 (CLM4.0), 0.72 ± 0.13 (NoahMP-BATS), and 0.68 ± 0.14 (NoahMP-CLASS) (Fig. 5a). Remotely sensed MODIS snow-covered albedo for cropland and grassland was lower (0.55 ± 0.11) than values averaged from the three land surface models.
The three land surface models overestimated snow-covered albedo of deciduous broadleaf forest relative to MODIS by 0.17 to 0.18 (Fig. 5a). Snow-covered albedo of mixed forest tended to be higher in CLM4.0 compared to the two Noah-MP options tested. For evergreen needleleaf forest, snow-covered albedo in NoahMP-BATS and NoahMP-CLASS was significantly lower than that from MODIS and CLM4.0 (Fig. 5a).
During snow-free dormant periods, all models agreed well with MODIS for deciduous broadleaf, mixed forest, and evergreen needleleaf forest (Fig. 5b). Cropland/grassland snow-free dormant albedo was higher in NoahMP-BATS and NoahMP-CLASS compared to MODIS and CLM4.0.
The relationship between modeled snow-covered albedo and snow-depth over cropland/grassland mosaic grid cells was compared to CoCoRAHS albedo and snow depth data collected over mowed lawns within the state of New Hampshire (Fig. 6). Modeled results were averaged by land surface model because results did not demonstrate any dependence on choice of microphysics or radiation schemes (not shown). For CLM4.0, the modeled results generally agreed well with the CoCoRAHS observations at snow depths greater than 5 cm (Fig. 6a). At snow depths less than 5 cm, CLM4.0 albedo was lower than CoCoRAHS measured albedo, likely due to the shorter canopy height of the CoCoRAHS mowed lawns (5–10 cm) compared to CLM4.0 short vegetation canopy height (20 cm). The snow albedo in NoahMP-BATS and NoahMP-CLASS agreed well with CoCoRAHS observations at snow depths greater than 10 cm and tends to be slightly higher in NoahMP-BATS (Fig. 6b) compared to NoahMP-CLASS (Fig. 6c).
b. Temperature responses to reforestation
1) Seasonal analysis of temperature responses
In all ensemble members, winter T2max warmed in response to replacing mid-1800s cropland/grassland with present-day forest cover. Ensemble members using the two Noah-MP land surface model options produced warming on the order of +3 to +10 K while the CLM4.0 simulations produced much weaker warming, ranging from +0.5 to +3 K (Figs. 7a,b). The warming response was stronger during the cold, snowy simulation (Figs. 7a and S5a) compared to the warm, dry simulation (Figs. 7b and S5b). The strongest warming in T2max occurred in the region of eastern, coastal Maine (i.e., “Downeast” Maine) where cropland/grassland was replaced with evergreen needleleaf forest (Figs. S5a,b).
The magnitude and sign of T2min responses to replacing cropland/grassland with present-day forest varied among ensemble members (Figs. 7c,d and S6). In general, replacing cropland/grassland with evergreen needleleaf forest resulted in a T2min cooling response that was slightly stronger in the warm, dry simulations than in the cold, snowy simulations (Figs. 7c,d and S6). Replacing mid-1800s cropland/grassland with deciduous broadleaf forest generally resulted in weak warming in T2min ranging from +0.2 to +1.5 K in the cold, snowy simulations (Fig. 7c) and from +0.2 to +0.6 K in the warm, dry simulations (Fig. 7d).
2) Diurnal analysis of surface energy fluxes
The differences in surface energy fluxes between the evergreen needleleaf forest and mid-1800s cropland/grassland were averaged by land surface model for each of the 3-h time steps during daylight hours. Differences in the daytime surface energy fluxes peaked around 15:00 local time (Figs. 8 and 9). For all three land surface models, replacing mid-1800s cropland/grassland with evergreen needleleaf forest resulted in an increase in shortwave radiation absorbed by vegetation (SABV), a decrease in shortwave radiation absorbed by bare ground (SABG), and an increase in sensible heat flux (HFX) (Fig. 8). Greater increases in SABV and HFX occurred in the cold, snowy simulations (Fig. 8a) compared to the warm, dry simulations (Fig. 8b). The increase in SABV was larger in the NoahMP-BATS and NoahMP-CLASS simulations compared to CLM4.0. Despite increases in sensible heat, the conversion of forest to cropland/grassland had little impact on planetary boundary layer height (PBLH; Figs. S7 and S8), convective available potential energy (CAPE; Figs. S9 and S10), or lifting condensation level (LCL; Figs. S11 and S12) in either the cold, snowy or warm, dry simulations.
At night, ground heat flux (GRDFLX) increased and HFX decreased when cropland/grassland was replaced with evergreen needleleaf forest (Fig. S13). Weaker responses in nighttime surface energy fluxes occurred in the cold, snowy simulations (Fig. S13a) compared to the warm, dry simulations (Fig. S13b).
When mid-1800s cropland/grassland was replaced with deciduous broadleaf forest, the differences in daytime (Fig. 9) and nighttime (Fig. S14) surface energy fluxes were smaller compared to replacement by evergreen needleleaf forest. During the day, the increase in SABV for deciduous broadleaf forest compared to cropland/grassland was 60–80 W m−2 (Fig. 9) compared to 220–300 W m−2 for evergreen needleleaf forest (Fig. 8). At night, GRDFLX decreased for deciduous broadleaf (Fig. S14), in contrast to evergreen needleleaf forest, which saw an increase in GRDFLX at night (Fig. S13).
This study used a multiphysics ensemble to evaluate the sensitivity of high-resolution regional historic land cover change simulations to land surface, radiation schemes, and microphysics schemes. The 12-member ensemble revealed that choice of land surface model and LWSW scheme produces substantially different biases in surface temperature, precipitation, and snow depth. For example, biases in T2max ranged from −4 K using WRF/CLM4.0 and CAM3.1 radiation to +8 K using WRF/Noah-MP and RRTMG radiation. For T2min, biases had little dependence on land surface model, but ranged from −6 K using CAM3.1 radiation to +4 K using RRTMG. The validation results indicated that no single “best” model configuration exists for simulations of winter climate in New England and, by extension, other mixed land-use temperate midlatitude regions. Use of a multiphysics ensemble provides a more informative and richer discussion of responses to historical land cover changes in the region than using a single model configuration.
a. Temperature sensitivity in the multiphysics ensemble
Our findings indicated a strong sensitivity of winter T2max to land surface scheme. A WRFV3.0 land surface scheme sensitivity study over the western United States from November 1995 through November 1996 also demonstrated that winter temperature was strongly influenced by choice of land surface model; however, sensitivity to radiation and microphysics schemes was not evaluated (Jin et al. (2010). Jin et al. (2010) compared four land surface schemes: 1) the simple soil diffusion scheme, 2) the original Noah scheme, 3) the Rapid Uptake Cycle scheme, and 4) CLM3.0. They found that CLM3.0 improved temperature simulations relative to other land surface schemes, with the exception of winter maximum temperature. In the current study, we find that CLM4.0 improves simulation of winter maximum temperature in the northeastern United States compared to NoahMP-BATS and NoahMP-CLASS. We note that the configurations of Noah-MP evaluated here are markedly improved relative to the original Noah scheme evaluated in Jin et al. (2010), as is CLM4.0 compared to CLM3.0. Nonetheless, WRF simulations exhibit a similar sensitivity to land surface model selection in both the current study and Jin et al. (2010).
A more recent WRFV3.5.1 sensitivity study evaluating the original Noah LSM, Noah-MP, and CLM4.0 over the western United States found that CLM4.0 outperformed Noah and Noah-MP in simulation of T2max, while Noah-MP performed best at simulating T2min (Chen et al. 2014). In our analysis, CLM4.0 was more skillful at simulating T2max than either of the two Noah-MP configurations we tested.
The current study demonstrated that winter T2min simulations were sensitive to choice of radiation scheme. Mooney et al. (2013) also observed sensitivity in winter (December–February) mean temperature to radiation scheme for simulations conducted over Europe, 1990–95; however, sensitivities to minimum temperature and maximum temperature were not evaluated separately.
We cannot draw strong conclusions about model performance between CLM4.0 and Noah-MP in the simulation of 2-m temperature given the differences in study regions, WRF Model versions, model resolutions, and physics schemes used in previous studies. However, we can conclude that a multiphysics ensemble helped identify model sensitivity to the choice of land surface scheme and radiation scheme in simulation of winter temperature.
b. Impacts of physics schemes on winter precipitation
Overall, the choice of land surface model and radiation schemes did not influence regionwide precipitation biases in the simulations. This is consistent with the findings reported in Jin et al. (2010), who compared four land surface models over the western United States but did not evaluate sensitivity of precipitation to other physics options (e.g., microphysics, radiation). Mooney et al. (2013) also identified a lack of sensitivity of European winter precipitation to land surface scheme and a general pattern of overestimation in areas of steep terrain that was also found in our study. Chen et al. (2014) found a 21%–26% overestimation of accumulated precipitation in western U.S. simulations and considerably better correlations (r2 ~ 0.92) with PRISM compared to the current study, where correlations (r2) were less than 0.6. Rasmussen et al. (2011) also found a general pattern of 10%–40% overestimation of November–May precipitation in Colorado WRF simulations relative to SNOTEL and PRISM, although they note that SNOTEL undercatch due to wind may reach up to 15%.
The overestimation of precipitation at high elevations in New England is challenging to diagnose given the dearth of ground-based observations at higher elevations in the region. Additionally, differences in study region, lateral boundary conditions, and physics options tested in previous studies make it difficult to explain why correlations between observed and WRF modeled precipitation over the New England region were considerably worse compared to the western U.S. WRF simulations. Only two configurations of Noah-MP were tested here; other combinations of Noah-MP options may result in improved simulation of winter precipitation.
c. Differences in albedo simulated by land surface models
Modeled snow-covered albedo over cropland/grassland (0.73 ± 0.15, 0.72 ± 0.13, and 0.68 ± 0.14 for CLM4.0, NoahMP-BATS, and NoahMP-CLASS, respectively) was generally much higher than MODIS albedo (0.55 ± 0.11). While MODIS averages for snow-covered cropland and grassland in New England agree well with global averages (0.58 ± 0.13 to 0.59 ± 0.11; Moody et al. 2007), MODIS albedo over cropland and grassland tends to be biased low compared to ground-based snow-covered observations over grassland at Fort Peck, Montana (0.75–0.90; Wang et al. 2014); cropland in Bondville, Illinois (0.7–0.85; Wang et al. 2014); pasture in Durham, New Hampshire (0.71; Burakowski et al. 2015); and grass lawns in New Hampshire (0.72; Burakowski et al. 2013).
The difference between modeled and MODIS albedo over grassland and cropland could be due to pixel heterogeneity within the 500-m resolution of MODIS data, whereas the WRF modeled albedo is reported for uniform land cover. Specifically, the IGBP land cover classification used here reports the predominant (>60%) land cover within a MODIS gridded 500-m pixel. Terrain undulations, roadways, buildings, vegetation protruding above the snowpack, or other contrastingly low albedo surfaces, could lower the snow-covered albedo within a 500-m MODIS pixel classified as grassland or cropland relative to smaller footprint of ground-based observations over relatively homogeneous surfaces (Liu et al. 2009).
Differences in the land models’ representations of snow-covered forest albedo could be due to either too much snow being represented on the canopy or too much snow visible underneath the canopy. In CLM4.0, canopy snow is an optical parameterization in which the snow-covered fraction of the canopy is used as a weight to average the scattering parameters used in the canopy with a temperature switch at 0°C (Oleson et al. 2010). The canopy snow temperature switch can lead to overestimation of snow-covered albedo over boreal forest canopies in CLM4.0 (Thackeray Fletcher and Derksen 2014), and could also impact deciduous broadleaf temperate forests. In Noah-MP, the snow-covered fraction of the canopy is similarly used as a weight to average the scattering parameters used in the canopy. However, the Noah-MP canopies have a maximum snow holding capacity that is dependent on vegetation-specific LAI, SAI, and the bulk density of the snowfall. In Noah-MP, canopy snow unloading responds to both vegetation temperature and wind (Yang et al. 2011). The low albedo bias of snow-covered evergreen needleleaf forest in Noah-MP relative to MODIS suggests that too much snow is unloaded from the canopy. This could have contributed to the warm biases in T2max in areas with evergreen needleleaf forest.
The comparison of modeled albedo and snow depth relationships revealed important differences in how the land surface models treat the burial of short vegetation. In CLM4.0, the LAI and SAI of short vegetation is adjusted until a critical snow depth of 0.2 m is reached, at which point the vegetation is considered completely buried by snow. In Noah-MP, the snow depth (hsnow,c) at which short vegetation (<0.5 m) is buried by snow is calculated as
where hυ,t is the canopy height and hsnow is the snow depth. When the modeled snow depth is greater than or equal to the critical snow depth, the fraction of buried vegetation is set to 1, and the effective LAI and SAI are set to zero. When the modeled snow depth is less than or equal to the critical snow burial depth, the fraction of buried vegetation is the snow depth divided by the critical snow depth and effective LAI and effective SAI are adjusted using the fraction of vegetation above the snowpack (Yang et al. 2011).
d. Climate responses to reforestation in the multiphysics ensemble
A few general patterns in T2max were observed when mid-1800s cropland/grassland was replaced with present-day forest: 1) evergreen needleleaf forest yielded the greatest warming response and deciduous broadleaf forest produced the weakest warming response, 2) the warming response was stronger in the cold, snowy winter simulations compared to the warm, dry simulations, and 3) the strength of the warming response in T2max was linked to the increase in shortwave energy absorbed by the forest vegetation and increase in sensible heat flux.
The multiphysics ensemble revealed a general agreement in sign, yet produced a wide range in the magnitude of responses to mid-1800s deforestation and subsequent regrowth of forests. At the high end of the spectrum, simulations using NoahMP-BATS and NoahMP-CLASS produced the strongest warming responses in T2max (+3 to +10 K). Warming simulated by Noah-MP ensemble members are unrealistic given an estimated warming trend of +1.5 K since 1850 in New England (Hodgkins et al. 2002), which includes the effects of land use change, anthropogenic emissions of greenhouse gases, and natural climate forcings (e.g., volcanic, solar). On the other end of the spectrum, simulations using CLM4.0 produced the weakest warming responses in T2max (+0.5 to +3 K), on par with the overall warming trend since 1850.
A number of factors contribute to the more realistic responses to reforestation in CLM4.0 compared to simulations using Noah-MP. Primarily, Noah-MP underestimated snow-covered albedo of evergreen needleleaf forests, which led to strong warm T2max biases relative to PRISM gridded observations. Second, both Noah-MP and CLM4.0 simulations had cool biases over open land relative to PRISM. When the warm-biased evergreen needleleaf forest in Noah-MP replaced the cool-biased cropland/grassland, the model response to reforestation in Noah-MP simulations was more pronounced than in simulations using CLM4.0.
The differences between replacing cropland with evergreen needleleaf versus deciduous broadleaf were broadly consistent with seasonal patterns identified using satellite-derived albedo, land surface temperature, and evapotranspiration by Zhao and Jackson (2014), who found that shortwave radiative forcing in winter was greater for the conversion of cropland to evergreen needleleaf than for conversion to deciduous broadleaf. This finding is consistent with the greater increase in shortwave radiation absorbed by vegetation reported here. Zhao and Jackson (2014) also derived radiative forcing due to changes in sensible heat and found that the conversion of cropland to evergreen needleleaf resulted in greater increases in sensible heat flux than the conversion to deciduous broadleaf, also consistent with the daytime results presented here (Figs. 8 and 9).
The warming in T2max when forest replaced cropland/grassland indicated that radiative forcing is the dominant biophysical effect during the day. Had energy redistribution from differences in surface roughness been the dominant factor, we would have expected the model to produce neutral responses in T2max when the forests replaced cropland/grassland. The neutral response results from the warming due to the radiative effect being offset by the cooling from the roughness effect, as forests can theoretically dissipate heat more efficiently than open land (Lee et al. 2011). We did not detect evidence of increased turbulence from forest canopies, as indicated by the insignificant changes in daytime PBLH, CAPE, and LCL when forest is replaced by cropland/grassland (Figs. S7–S12). Simulation of finescale turbulence remains challenging but would help improve understanding of how land use affects surface climate.
For T2min, the range in the multiphysics ensemble response was much smaller compared to the range in responses in T2max and the sign of the temperature response varied depending on what type of forest replaced the mid-1800s cropland/grassland. All ensemble members produced warming at night in T2min when present-day deciduous broadleaf replaced cropland/grassland, a result that was associated with a decrease in ground heat flux. The increase in ground radiative heat flux detected when evergreen needleleaf replaced cropland/grassland indicates that more heat was lost from the ground, leading to cooling at the surface. Ground radiative heat flux is generally negative at night (heat is lost from ground to the land surface–atmosphere interface), and more negative for grassland compared to forest at night (Oliver et al. 1987). When snow cover is present it insulates the ground surface and reduces heat loss from the ground.
The warmer T2min (+0.5 to +1.5 K) over forest compared to cropland/grassland reported in our study is consistent with tower observations presented by Lee et al. (2011), who reported that temperate forests were about 1.5 K warmer than adjacent open lands at night, primarily due to energy redistribution from changes in surface roughness. However, we identified cooling in T2min (0.5 to 2.2 K) over evergreen needleleaf forest compared with cropland/grassland that was in contrast to Lee et al.’s (2011) tower observations. It is worth noting that the temperature responses to deciduous broadleaf and evergreen needleleaf in the model simulations presented here were driven by changes in nighttime ground heat flux; there was no impact on nighttime turbulence that could be detected by comparing differences in PBLH, CAPE, and LCL (Figs. S7–S12). This indicated either that the energy redistribution from changes in surface roughness and partitioning of sensible and latent heat were not well simulated by the model or that more observations are needed to widely confirm nocturnal warming from reforestation. We also note that the multiphysics ensemble presented here did not include ensemble testing of turbulence schemes or PBL schemes. Future work should evaluate sensitivity to these schemes and further investigate the effects of energy redistribution from changes in surface roughness.
The winter climate impacts of large-scale deforestation and subsequent reforestation in the New England region of the United States and other midlatitude temperate regions have received relatively little attention. Here, we used a multiphysics ensemble of regional climate simulations to investigate responses of surface climate to cropland abandonment and subsequent reforestation.
The multiphysics ensemble revealed sensitivity in T2max to choice of land surface scheme and sensitivity to radiation scheme in T2min. The largest T2max warm biases occurred in ensemble members using Noah-MP, specifically in regions with low biases in snow-covered evergreen needleleaf albedo relative to MODIS and CLM4.0. None of the ensemble members simulated regional precipitation adequately enough to evaluate the impact of land cover change on precipitation patterns.
The climate response to mid-1800s deforestation and subsequent reforestation in New England very likely led to warming in daytime T2max. Much of the warming in T2max can be attributed to a decrease in albedo and subsequent increase in shortwave radiation absorbed by the forest canopy compared to cropland/grassland. The decrease in surface albedo increased the sensible heat flux at the surface, leading to warmer temperatures during the daytime. The magnitude of T2max warming remains uncertain. Because of the large warm biases in T2max in the Noah-MP ensemble members, the range in T2max warming (+0.5 to +3 K) in the CLM4.0 simulations is likely a better approximation of the winter climate response to reforestation that occurred in the New England region from the mid-1800s to present day. The weaker T2max warming responses in the warm, dry simulations suggest that future projections of warmer and lower snowfall winters could diminish radiative forcing associated with future changes in forest cover.
The T2min response at night remains more uncertain; only a few of the multiphysics ensemble members captured the nocturnal warming identified in surface observations (Lee et al. 2011). The ensemble members that did capture warming did so only over deciduous broadleaf forest and some areas of mixed forest. Over evergreen needleleaf forest, surface temperatures cooled relative to open lands. Future research should evaluate in greater detail the energy redistribution associated with surface roughness over a variety of forest types, as this appears to contribute significantly at night in observations and its magnitude was not currently captured in the multiphysics ensemble evaluated here.
The authors acknowledge the support of the National Science Foundation Experimental Program to Stimulate Competitive Research (EPSCoR), the New Hampshire EPSCoR Ecosystems and Society Program (EPS-1101245), the NASA Terrestrial Ecology Program (NNX12AK56G S01), the Northern Ecosystems Research Cooperative (NERC), and National Science Foundation Grant EF-1048481. This study is contribution of the NSF-supported Harvard Forest and Hubbard Brook Long-Term Ecological Research projects. We would also like to recognize the computational support provided by the NCAR Computational and Information Systems Laboratory for computing, processing, and data storage. We thank two anonymous reviewers for their contributions toward improving the manuscript.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/JCLI-D-15-0286.s1.
The National Center for Atmospheric Research is sponsored by the National Science Foundation.