Linking global land surface temperature projections to radiative effects of hydrometeors under a global warming scenario

Land skin temperature (Ts) is directly influenced by surface energy balance, in particular, radiative energy, which can be linked to model’s representation of radiative effects of hydrometeors in the atmosphere. This link is inferred by examining the changes of geographical distribution and seasonal cycle of surface radiation, surface turbulent fluxes and Ts between a pair of 140 years sensitivity experiments under 1% per year increase of atmospheric CO2. One is with radiative effects of falling ice (snow) hydrometeors on (SON) and the other off (NOS) using CESM1-CAM5 and the results are compared with CMIP5 models without these effects. For boreal winter, NOS relative to SON simulates less surface downward longwave and net flux (∼10–15 W m−2), resulting in colder Ts (∼2–3 K colder), over mid- and high latitudes, but more solar radiative flux, resulting in warmer Ts (∼1–3 K), over subtropical and tropical land. These differences between NOS and SON are amplified as the surface and the atmosphere become warmer. The results from CMIP5 ensemble generally match with those of NOS. Temporal correlation analysis indicates that the linkage between Ts and falling ice hydrometeor changes is through one between Ts and downward longwave and net fluxes at high latitudes, but strongly weakened by shortwave changes at low latitudes (and boreal summer). Relative to SON, land skin temperatures in NOS and CMIP5 are underestimated throughout the seasonal cycle but only slightly in summer.


Introduction
Land surface energy balance and its changes are parts of a complex land-atmosphere interactive system. They are a key factor for the climate change, which is influenced by the net surface radiative fluxes at a variety of temporal and spatial scales (Susskind 1993, Jin et al 1997, Pitman 2003, Knutson et al 2013, Brooks et al 2015. The land radiative skin (surface) temperature is primarily controlled by the atmospheric radiative fluxes at the skin interface, which acts as an upper boundary that constrains the land surface energy and hydrological cycle and then influences surface turbulent fluxes (latent and sensible) and surface precipitation (Gu et al 2006, Seneviratne et al 2006, 2010, Koster et al 2009, Schwingshackl et al 2017, Haghighi et al 2018.
Despite that land-related processes within land surface models have been advanced with the help of satellite measurements of, for example, albedo, emissivity, and leaf area index (Pitman 2003), they still suffer from the uncertainties in representations of physical processes from the atmospheric components of coupled global climate models (GCMs) such as clouds, precipitation and their interactions with radiation (Stephens 2005, Randall et al 2007, Trenberth and Fasullo 2009, Stephens et al 2012a, Stephens and Li et al 2012b. Land surface properties change with those of the upper boundary conditions from atmospheric radiation forcing.
It is therefore critical to investigate the impacts from the upper boundary conditions related to clouds and radiation on land surface properties (Lo and Famiglietti 2013, Anderson et al 2015, Wey et al 2015. As discussed below, the lack of the fidelity in the presentday CMIP (Coupled Model Intercomparison Project) model land surface simulations can lower the confidence level for future projection of land surface properties (Hua et al 2014, Li et al 2016a, 2016b. It has been reported that modeled land surface temperature of the present-day climate simulations in CMIP Phase 5 (CMIP5; Taylor et al 2012) seems to be underestimated in mid-and high-latitudes in boreal winter season but biased warm at low-and midlatitudes, especially over arid and semiarid regions in boreal summer season (e.g. Morcrette 2018), attributed mainly to the discrepancies of the net shortwave radiation and surface longwave radiation (Hua et al 2014, Li et al 2016a, 2016b. That is, overestimation of downward shortwave radiation at the surface (RSDS) warms the land surface, mostly during boreal summer, while underestimation of downward longwave radiation (RLDS) cools the land surface over most continental regions, in particular, during boreal winter. Li et al (2016a) discussed the key role of treating radiative effects of falling ice (snow) hydrometeors in simulating the present-day seasonal cycle of land skin temperature variations. They found that such a treatment is more important for the boreal winter season than the other seasons. This is mainly because downward longwave radiation increases due to the greenhouse effect of falling ice hydrometers, which in turn restricts the cooling of land surface when there is less solar radiation during winter.
Since land surface temperature and surface energy budgets are critical factors that control near surface air temperature through surface layer turbulent fluxes, their biases in most CMIP5 GCMs, which treat the radiative effects of floating ice hydrometeors but not those of falling ice (snow) hydrometeors, impose non-trivial impacts on modeled land-atmosphere interactions, soil layer and land hydrological processes through soil moisture and temperature changes (Li et al 2016b). These biases in the present-day climate simulations may influence the confidence level for projecting future land surface properties. Only a few CMIP5 GCMs treat the radiative effects of both floating ice and precipitating ice hydrometeors, including the National Center for Atmospheric Research-Department of Energy (NCAR-DOE) CESM1-CAM5 (Neale 2012) and HadGEM2 (Martin et al 2011). In this study, we choose the former, which is a community model, for a sensitivity study.
Following Li et al (2016a), in this study, we will use CMIP5 outputs to analyze the linkage of the projected land surface temperature with radiative effects of hydrometeors via surface energy budget. We focus on the physical process related to the falling ice (snow) radiative effects (FIREs) using controlled CESM1-CAM5 simulations, with an emphasis on the geographic distribution changes of projected land surface temperature and the associated relationships with surface energy budget components under a global warming scenario. We discuss how inclusion of FIREs, compared to exclusion of FIREs, can substantially change the simulated surface energy budgets and land surface temperature focusing on boreal winter (December, January and February, DJF), which is chosen due to the strongest signal in the seasonal cycle, and how the seasonal cycle of changes is impacted by FIREs for a few selected regions.
In section 2, we describe the controlled simulations and analysis methods. Results are presented in section 3. We discuss and conclude major findings in section 4.

Controlled climate model simulations
Two simulations are performed, one with FIREs (snow ON or SON), in which radiative effects of both floating ice and falling ice hydrometeors are included, and another without FIREs (NO snow or NOS), in which radiative effects of floating ice hydrometeor only are included, using CESM1-CAM5 (Neale 2012). Following the CMIP5 protocol, both climate change simulations are initialized from a preindustrial control (piControl) with 1% per year increase in atmospheric CO 2 , resulting in a quadrupling of CO 2 at the end of 140 years simulation, to examine the responses of land surface under a gradually warming climate. Radiative effects of cloud liquid hydrometeor are included in both. Similar simulations with 13 CMIP5 models without FIREs are considered (table 1), and their ensemble mean is abbreviated as CM5. 'CM5-SON' and 'NOS-SON' will be used in the rest of this paper to denote the differences of CM5 and NOS from SON, respectively.
A detailed documentation of the Community Earth System Model Version 1 (CESM1) is available at www.cesm.ucar.edu/models/cesm1.0/, which is a fully-coupled earth system model that consists of the Earth's atmosphere, ocean, land, sea ice and land ice (glacier) components. We only provide a brief description of the stratiform cloud microphysics parameterization used in the atmosphere component (CAM5: Community Atmosphere Model Version 5) here. This parameterization is based on the Morrison-Gettelman scheme (MG1: Morrison and Gettelman (2008)), which is a four-water species (liquid, ice, rain, and snow), two-moment cloud microphysics. Cloud water and ice are prognostic but rain and snow are diagnostic. The falling ice large particles (also known as snow) are represented in the model with appreciable falling velocities and radiatively active using the diagnosed mass and effective radius Gettelman 2008, Gettelman et al 2010). The simulated cloud ice (also known as floating ice) and falling ice (accounting for about two third of total ice water path) are comparable to the CloudSat retrieval (Gettelman et al 2010). Further details regarding the diagnosed snow mass is described in Morrison and Gettelman (2008) and its radiative properties in CAM5 are provided in Gettelman et al (2010).

Analysis method
The surface energy balance and the relationships of its components with land surface temperature are analyzed following the method presented in an earlier study (Li et al 2016a). The energy balance components include surface downward shortwave flux (RSDS), surface downward longwave flux (RLDS), reflected shortwave flux (RSUS), latent heat flux (LHFX) and sensible heat flux (SHFX), as well as ground net heat flux (GHFX). A small and slowvarying net ground heat uptake or release from the soil layers is a good approximation for studying nontransient behavior (20 years average, in this study). With this assumption, the skin radiative (or land surface) temperature (Ts) is directly linked to the modeled land surface energy budget components, based upon the surface energy balance. We briefly describe the land surface energy balance below. More details are referred to Li et al (2016a). The local land surface temperature (Ts) can be derived from the Stefan-Boltzman law shown in equation (1): where Femit is the longwave radiation emitted from the surface, ε is the grid-box averaged bulk emissivity, and σ is the Stefan-Boltzman constant given by The land surface temperature (Ts) can be linked to surface energy balance by equation (2): where a is the grid-box averaged albedo, and upward flux from surface to the atmosphere (SHFX, LHFX, GHFX) is defined as positive. GHFX is generally very small in the long-term time average over a GCM gridbox: Assuming that emissivity, ε, is close to unity and RSUS = a RSDS, equation (3) can be approximately expressed as: The combination of the three radiative fluxes with surface turbulent heat fluxes on the right-hand-side of (4) represent the response of surface thermal emission to the net heat absorbed at the land surface, that is, RLUS = σT 4 s , which can be called 'net radiative flux (Net).' All the fluxes and Ts are provided by the climate models, which allow us to explore the direct linkage of Ts with surface radiation fluxes and surface latent and sensible heat fluxes, and indirect linkage of Ts with radiative effects of hydrometeors by comparing the results between NOS and SON simulations.

Changes of boreal winter hydrometeor water paths
To understand the linkage of surface energy budget and Ts to radiative effects of hydrometeors, we first examine three hydrometeor water paths (liquid water path or LWP, ice water path or IWP, and snow water path or SWP) and their changes that directly impact the radiative energy changes (figure 1). We first show the geographical distribution of the difference (NOS-SON) between NOS and SON in search of direct and indirect linkages of LWP/IWP/SWP changes with changes in radiation (RLDS and net flux) and Ts. Because there is no impact of SWP on radiation in NOS, which is thus assigned to be zero, negative values of SWP and its changes between the first (1-20) years and last (121-140) 20 years in SON are shown in figures 1(a)-(c).
The geographical distributions of SWP in SON are similar to the present-day observations (Li et al 2012) for the initial 20 years mean state (first period), with large values (60-150 g m −2 ) over Europe, North America as well as south Africa and Central America ( figure 1(a)). Relative to the first period, the spatial pattern of SWP for the last 20 yr mean state (last period) remains approximately the same but with larger magnitudes in some regions ( figure 1(b)). The changes in SWP from the first to the last periods occur over midlatitudes and high latitudes of Eurasia and North America, with a magnitude of about −5 to −30 g m −2 (figure 1(c)), meaning that SWP in SON increases over these latitudes as the surface and the atmosphere become warmer. Small areas with an opposite sign of changes occur over the lower latitudes (except for central Africa). Except for LWP changes over Europe (10-15 g m −2 ), changes in LWP and IWP (2-10 g m −2 ) over small areas between the two periods from NOS to SON (figures 1(f) and (j)) are opposite in sign to those in SWP, meaning that changes in SWP mostly neutralize those in LWP/IWP. This result clearly suggests that changes in hydrometeors between the two periods are most significant over the midlatitudes and high latitudes in boreal winter, especially for SWP, which may imply more impacts on radiation and Ts from falling ice in the SON simulation than from floating ice or liquid water (except for a small area over Europe).

Changes of boreal winter surface energy fluxes and land surface temperature
In this section, we search for direct linkages between surface energy fluxes and Ts and indirect linkage between SWP and Ts from the geographical distributions of NOS-SON differences for the first and last periods and their changes between the two periods.
The following question will be addressed: does the exclusion of FIREs in NOS relative to SON cause less RLDS, more RSDS and RSUS, and thus less net flux, resulting in colder Ts under the progressive global warming, as documented in Li et al (2016a) for the present-day climate? Figure 2 shows the NOS-SON differences in surface energy fluxes, Ts and the ratio of Ts to net flux for the first period. The geolocations of maximum SWP differences (figure 1(a)) are in general agreement with less net flux (figure 2(a)), which are mainly contributed by less RLDS (figure 2(b)), and with colder Ts (figure 2(d)) over the entire mid-and highlatitude north of ∼40 • N, in particular, over the western part of North America and western Europe and Antarctic continent. More RSDS (figure 2(f)) with 5-10 W m −2 between the latitude belts of 60 • N-35 • S, where the solar flux dominates (figure 2(a)) such as over North America and Asia, South America, Australia and Africa with warmer Ts by 1-2 K south of 25 • N (figure 2(d)). Latent heat flux (LHFX: figure 2(h)), sensible heat flux (SHFX: figure 2(e)) and RSUS (figure 2(i)) are generally rather small (<5 W m −2 ) globally compared to other radiative fluxes. However, there are appreciable differences in SHFX over the regions north of 35 • N, in LHFX over the Southern Hemisphere, and in RSUS over northern midlatitudes and the Antarctic. The mean Ts per unit net flux is not small for the regions north of 35 • N, Greenland and the Antarctic continent, about 0.1-1 K per W m −2 , as shown in figure 2(g), suggesting a strong direct linkage between Ts and net flux over many land areas.
The geolocation distributions of NOS-SON differences in surface energy fluxes and Ts for the last period are similar but with larger magnitudes of differences (figure A1), shifting towards the higher northern latitudes, and broader spatial extents compared to those for the first period (figure 2). The NOS-SON Ts differences become more negative, i.e. much colder, in winter from the first to last periods (figure A1(d)) with a magnitude of 2-3 K, which is directly related to less RLDS (figure A1(b)) and net flux ( figure A1(a)).
As in figures 2 and A1, the geographical patterns of changes in the NOS-SON differences from the first to last periods are shown in figure 3. At first glance, the geolocation distributions are rather similar to those in the first (figure 2) and last periods (figure A1), except for those of the mean Ts per unit net flux (figure 3(g)) and the latitudinal shifting of the patterns. However, changes in Ts are more negatives at higher northern latitudes (∼40 • N), slightly more positive between the equator and 40 • N ( figure 3(d)). This means that NOS relative to SON tends to produce colder Ts up to 2-3 K over the high latitudes as the surface and the atmosphere become warmer, which is directly related to the NOS-SON changes towards less net and RLDS fluxes (figures 3(a) and (b)), especially over western Europe and northern and eastern America.
Over the regions south of ∼40 • N including the entire Southern Hemisphere, however, the NOS-SON changes in Ts are either very small or slightly positive. This is due to the combination of positive changes in RLDS and RSNS (net solar radiation; figures 3(b) and (c)) from the first to last periods. The negative changes in LHFX in the Southern  Hemisphere (figure 3(h)) and those of SHFX over some small areas (figure 3(e)) also contribute. The changes in mean Ts per unit net flux (figure 3(g)) are relatively small north of 40 • N. Elsewhere, there are large positive and negative changes over small areas, which is contributed by downward longwave (due to increasing LWP) in boreal winter (figure 3(b)) but downward shortwave (RSDS) in austral summer ( figure 3(f)).
The strong similarity in geolocations between Ts and surface energy flux changes reveals their direct linkages for the changes between the two periods. To further understand the linkages, we detrend NOS-SON Ts, RLDS, RSDS, RSUS and net flux and show the geographical distributions of the temporal correlation of radiative flux with Ts differences for DJF using the 140 years output in figure 4. The NOS-SON differences in Ts have very high positive correlations with those of RLDS over most of Europe, North and South America, Africa and Australia, but negative correlations with RSDS and RSUS over the northern mid-and high-latitudes and the Antarctic and weak or positive correlations over the low latitudes and Southern Hemisphere continents. The net solar flux (not shown) is positively correlated with Ts over the low latitudes and Southern Hemisphere continents, which leads to the positive correlations of Ts with net flux everywhere. Further, SWP is positively correlated with RLDS and Ts over the northern mid-and high-latitudes and the Antarctic, but negatively correlated with RSDS and RSUS over low latitudes ( figure A2). This result may suggest that the indirect linkage between Ts and SWP is through the linkage between Ts and RLDS at high latitudes (boreal winter), but strongly weakened by shortwave changes at low latitudes (and boreal summer) (figure 4(b)) and due possibly to changes in other hydrometeors (figures 1(f) and (i)).
In summary, relative to NOS, SON produces consistently higher net downward surface radiative fluxes (figure 2(a)) and its changes (figure 3(a)) due to increased downward longwave heating contributed by increased SWP in winter, which increases Ts. In summer, changes in net radiative flux can be either positive or negative, as implied from the results of the Southern Hemisphere lands, depending on the offset between increased RLDS and decreased RSDS due to FIREs, which influences the changes in Ts (not shown). This is also true for the lower latitudes in winter due to the increasing role of solar radiation. These results suggest that there is an indirect linkage between changes in Ts and SWP over high latitudes. The linkage is weak over low latitudes due possibly to contributions of an opposite sign from other hydrometeors (figures 1(f) and (i)).

Comparisons of CMIP5 models with CESM1-CAM5
Based upon the credential of using NCAR CESM1-CAM5 with SON relative to NOS in simulating the present-day climate, which shows an improved simulation of radiation fields and Ts (Li et al 2016a), we use SON as a reference for comparisons with ensemble model mean of CMIP5 models (table 1), i.e. CM5. Similar to figure 3 for NOS-SON, figure 5 shows geographic distributions of the changes of CM5-SON differences between the first and last periods. SHFX and LHFX changes are small and thus not shown. The patterns of changes in RLDS, RSDS, RSUS, Net, and Ts are similar to those in NOS-SON (figure 3) except over the subtropical and tropical land areas due to the CM5-NOS differences (figure A3). It is important to point out that CM5 projects larger negative Ts changes (magnitudes of −2 to −3 K) over Europe, north of 40 • N between 30 • E and 90 • E and over Alaska and west Canada, which is reinforced by the same sign of the CM5-NOS differences (figure A3), but not the eastern part of North America (i.e. warming). The changes in each of the four radiative fluxes are generally similar between NOS-SON and CM5-SON but with reduced impacts of FIREs over the high latitudes and enhanced CM5-SON differences in regions with small NOS-SON differences, due to the CM5-NOS differences ( figure A3). This result further illustrates that RLDS is important in determining the relationship between net flux and Ts for DJF ( figure 5(g)), from which we can further infer an indirect linkage between changes in SWP and Ts. The above broader regions disparities seem to be partially contributed by the basic differences of CESM1-CAM5 relative to CMIP5 models, which can be resulted from the large model spread in simulating the regional characteristics of progressive climate changes among CMIP5 models (Li et al 2016a), e.g. cooling extends from North to Central America ( figure 5(d)). In addition, differences in model physical processes other than FIREs may also contribute to the large model spread (as shown below) and different projections in different regions.

Long-term trends and changes of the seasonal cycles
Shown in figure 6 is the time series of boreal winter (DJF) Ts and radiation components for NOS, SON and CM5 averaged over the 40 • N-70 • N latitude belt but excluding the Greenland where land ice dominates. SON shows higher Ts by ∼1-2 K (figures 6(a) and (b)), accompanied by larger RLDS (6c) and Net flux (6f) but smaller differences in RSDS (6d) and RSUS (6e), compared to NOS and CM5 for the entire 140 years. Although NOS and CM5 are largely similar to each other, figures A6(a) and (b) indicate that CM5 has a slower warming rate than either NOS or NOS especially after year 70, due to the differences of NOS relative to CM5, suggesting that besides the inclusion of FIREs, other physical processes in CMIP5 models are likely involved. It is noted that Ts trends are nearly identical and so are those of RLDS and Net flux between NOS and CM5 before year 50. Over Europe (50 • N-80 • N, 50 • E-70 • E; figure A4) and North America (140 • E-80 • W, 40 • N-70 • N; figure A5), the results are similar to those shown in figure 6 except that the differences of NOS relative to CM5 are larger over Europe, as noted in figure A3.
We further compare seasonal cycle of Ts and radiation components for SON against NOS (NOS-SON) and CM5 (CM5-SON). For this, we select the Eurasian region (0-180 • E, 40 • N-70 • N), where Ts and radiative energy flux changes are strongly correlated (figures 3-5). The annual cycle of ensemble mean and standard deviation are calculated from each CMIP5 model with its annual-mean being removed. The annual cycles of CM5 model ensemble mean (thick black) with one plus/minus standard deviation (thin black) relative to SON are plotted in figures 7(a)-(e), along with NOS-SON (red). The time variations of each parameter in the first (figures A6(a)-(e)) and last periods (figures 7(a)-(e)) are similar to each other with colder Ts, less RLDS and Net in NOS relative to SON. Thus, we only present the results for the last period and the differences between the last and first periods. Figures 7(a)-(e) show that NOS-SON Ts is colder by up to −3 K during September-April, along with negative RLDS and Net flux differences up to −10 W m −2 , but it is slightly warmer during spring-summer, accompanied by positive net flux differences due to RSDS ( figure 7(c)). Although the Ts amplitude of CM5-SON seasonal cycle is small, the amplitudes of radiative energy fluxes are much larger (up to 30 W m −2 ), compared to NOS-SON, along with large standard deviation (up to 4 K for Ts, 20 W m −2 for radiative fluxes), in particular, in RSDS (figure 7(c)) and net flux (figure 7(e)), indicating large inter-model differences. There are also phase differences in the seasonal cycles between CM5-SON and NOS-SON. These differences between SON-NOS and CM5-SON may be attributed to physical processes and inter-model spread other than FIREs. For example, variations of SHFX, LHFX, ground wetness and evaporation and vegetation evapotranspiration may adjust to large surface radiative energy imbalances differently in different models so that the seasonal cycle of Net flux does not solely determine that of Ts.
Figures 7(f)-(j) show that both NOS and CM5 project colder Ts (up to −2 K) than SON throughout the entire seasonal cycle, with phase differences between NOS-SON and CM5-SON. The changes in RLDS and net flux (figures 7(g) and (j)) are mostly consistent with those in Ts except for larger amplitudes in CM5-SON, due to positive changes in RLDS for October-November and large positive RSDS between January and July (figures 7(h) and (i)). The much large net flux changes in CM5-SON may suggest that large inter-model spread and SHFX and LHFX may play a more important role in regulating the annual cycle of individual CM5 models.

Discussion and conclusions
The linkage of global climate warming projections over land surface with radiative effects of hydrometeors has been explored by examining geographical distribution changes of radiation, turbulent fluxes and land surface temperature (Ts) and their seasonal cycles under the warming scenario of 1% CO 2 increase per year for 140 years. Two simulations are performed, one with FIREs (snow ON or SON), in which radiative effects of both floating ice and falling ice hydrometeors are included, and another without FIREs (NO snow or NOS), in which radiative effects of floating ice hydrometeor only are included, using CESM1-CAM5. Radiative effects of cloud liquid hydrometeor are included in both. Similar simulations with 13 CMIP5 models without FIREs are also compared with SON and NOS. As in Li et al (2016a), the surface energy balance is used to examine the direct linkage of Ts with surface radiation fluxes and surface latent and sensible heat fluxes, and infer an indirect linkage of Ts with radiative effects of hydrometeors by comparing the results between NOS and SON simulations.
For boreal winter, NOS, relative to SON, simulates less surface downward longwave and net flux (∼10-15 W m −2 ), resulting in colder Ts over midand high-latitudes (∼2-3 K colder), in particular, western Europe and northeast America, but more solar radiative flux, resulting in warmer Ts over subtropical and tropical lands (∼1-3 K). These NOS-SON differences are amplified in some regions as the surface and the atmosphere become warmer. We then examined geographical distributions of changes in terms of boreal winter (DJF) means and seasonal cycle of radiative fluxes (RLDS, RSDS, RSUS, and net flux) and Ts from the first (years 1-20) to last (years 121-140) periods. The changes of NOS-SON differences in Ts are more negative at higher northern latitudes (∼40 • N), slightly more positive between the equator and 40 • N. This is directly related to the NOS-SON changes towards less net and RLDS fluxes over western Europe and northeast America. Over the regions south of ∼40 • N including the entire Southern Hemisphere, however, the NOS-SON changes in Ts are either very small or slightly positive. This is due to a combination of positive changes in RLDS and RSNS (net solar radiation) although sensible and latent heat fluxes also matter. The strong similarity in geolocations between Ts and surface energy flux changes reveals their direct linkages for the changes between the two periods, as also confirmed by the temporal correlation analysis and significance tests (figure A7) using the 140 yr data. The results also suggest a stronger indirect linkage between Ts and SWP changes over the high latitudes than over the low latitudes.
During summer (JJA), the situation is more complex (figures not shown; but similar to boreal winter Southern Hemisphere), due to the cancellation of FIRE impacts between increased RLDS and decreased downward SW, resulting in non-systematic changes of Ts cooling/warming over small areas. The annual mean change (figure A8) is similar to that in winter (DJF) except for smaller magnitudes. The increased falling ice mass in SON leads to increased RLDS, resulting in a higher warming rate in Ts during the 140 year integration, compared to the NOS simulation.
We also compared the NOS-SON changes in Ts and radiative fluxes with those of CM5-SON (CM5: CMIP5 models) changes between the first and last periods, which are very similar except over the subtropical and tropical land areas. Both NOS and CM5 project higher changes in cooling (magnitudes of ∼2-3 K) over broad areas of Europe, north of 40 • N between 30 • E and 90 • E and over Alaska and west Canada, but not the eastern part of North America (i.e. warming), due to the basic differences between NOS and CM5, which seem to be the product of ensemble average resulted from the large model spread in simulating the regional characteristics of progressive climate changes among CMIP5 models (Li et al 2016a). Differences in model physical processes other than FIREs may also contribute to the large model spread as the response of surface sensible and latent heat fluxes to radiative energy imbalance may be different among models. Further analysis of the cause of these differences is beyond of the scope of this study.
As discussed above, geographic patterns of change are mostly similar between CM5 and NOS and so are the seasonal cycles. This resemblance implies that neglecting the falling ice-radiation interaction contributes, at least partially, to too cold mid-and high latitudes Ts under the CO 2 -driven global warming. We found that with FIREs, the contributions of net surface flux are primary from surface downward longwave flux for all seasons. The associated surface downward shortwave flux and surface reflected shortwave flux also play a non-negligible role during the summer season but somewhat offset by downward longwave counterpart. In these CESM1-CAM5 simulations, it appears that the inclusion of FIREs drives faster land surface warming over mid-and highlatitudes due to the fact that the greenhouse effects of snow hydrometeors affect closer to the land surface than over low latitudes.
We acknowledge the extensive studies from many recent studies (e.g. Singh et al 2019, Lehner et al 2020) regarding the potential contribution from internal climate variability, which could be one of the factors that influence the potential impact of the FIREs on land surface change projection presented in this study. While we have found the significant impacts from the FIREs on the changes of projected land surface radiation budgets, skin temperature, there are caveats and limitations for not considering the internal variability's contribution to the climate changes due to using a single model simulation in this study (Lehner et al 2020).
The results presented in this study have shown that radiative effects of snow are directly linked to the changes of surface energy balance and indirectly linked to Ts over the land regions under 1% CO 2 increase per year scenario. Other physical-dynamical processes have not been examined in this study and may also play a substantial role in Ts projection. Nevertheless, this indirect linkage is missing in almost all CMIP5 models and a majority of CMIP Phase 6 models (Eyring et al 2016), which project lower land surface warming over the high latitudes as in NOS. If these radiative effects are included across all models, it would potentially increase our confidence on the projections of regional climate changes and climate feedbacks. Inferring from results found in this study, it is difficult to conclude that the inclusion of FIREs is the only reason for an improved projection unless tests similar to CESM1 NOS and SON are performed with other models to isolate the FIREs, which is beyond the scope of this study.

Data availability statement
No new data were created or analyzed in this study.