Deciphering the sensitivity of urban canopy air temperature to anthropogenic heat flux with a forcing-feedback framework

The sensitivity of urban canopy air temperature ( Ta ) to anthropogenic heat flux ( QAH ) is known to vary with space and time, but the key factors controlling such spatiotemporal variabilities remain elusive. To quantify the contributions of different physical processes to the magnitude and variability of ΔTa/ΔQAH (where Δ represents a change), we develop a forcing-feedback framework based on the energy budget of air within the urban canopy layer and apply it to diagnosing ΔTa/ΔQAH simulated by the Community Land Model Urban over the contiguous United States (CONUS). In summer, the median ΔTa/ΔQAH is around 0.01 K W m−2−1 over the CONUS. Besides the direct effect of QAH on Ta , there are important feedbacks through changes in the surface temperature, the atmosphere–canopy air heat conductance ( ca ), and the surface–canopy air heat conductance. The positive and negative feedbacks nearly cancel each other out and ΔTa/ΔQAH is mostly controlled by the direct effect in summer. In winter, ΔTa/ΔQAH becomes stronger, with the median value increased by about 20% due to weakened negative feedback associated with ca . The spatial and temporal (both seasonal and diurnal) variability of ΔTa/ΔQAH as well as the nonlinear response of ΔTa  to ΔQAH  are strongly related to the variability of ca , highlighting the importance of correctly parameterizing convective heat transfer in urban canopy models.


Introduction
Anthropogenic heating resulting from energy consumption by human activities is an important controllor of urban climate. Although they occupy only 3% of the Earth's surface, cities consume 60%-80% of global energy and house more than half of the human population (United Nations 2022). The intense anthropogenic heating in cities can increase heat stress (Doan et al 2019, Jin et al 2020, Molnár et al 2020, which threatens thermal comfort and causes heat-related illnesses (Mora et al 2017). Studies have found that air warming of 1 • C is associated with a 1.8% increase in the mortality rate in cities when the daily temperature is higher than 28 • C (Chan et al 2012). Meanwhile, a higher temperature resulting from anthropogenic heating affects cooling energy demand, air quality, ecosystems and so on (Fink et al 2014, Xie et al 2016, Liu et al 2020. Anthropogenic heat flux also affects meteorological processes within the urban boundary layer (Fan and Sailor 2005, Chen et al 2009, Suga et al 2009, Krpo et al 2010, Bohnenstengel et al 2014, Zhang et al 2016, Ma et al 2017, Molnár et al 2020, Mei and Yuan 2021. Anthropogenic heat flux is generated from many sources, including building and industrial energy consumption, traffic and human metabolism (Sailor 2011, Chow et al 2014, Sun et al 2018. The magnitude of anthropogenic heat flux varies strongly with the local climate, population density, economy and technology (Fan and Sailor 2005, Allen et al 2011, Sailor et al 2015, Yang et al 2017, Jin et al 2020. The magnitude of anthropogenic heat flux is also scale dependent. At long-term and city (or larger) scales, the anthropogenic heat flux is typically of the order of 0.1-1 W m −2 . For example, Sailor et al (2015) developed a national database of anthropogenic heat flux over the contiguous United States (CONUS) and showed that the maximum wintertime (summertime) anthropogenic heat flux is around 0.8-0.97 W m −2 (0.47-0.63 W m −2 ) across 61 US cities. Another study reported that the annual mean anthropogenic heat flux is around 0.39 W m −2 , 0.68 W m −2 and 0.22 W m −2 for COUNS, western Europe and China, respectively, and only 0.028 W m −2 on the global scale (Flanner 2009). However, the short-term and neighborhood-scale anthropogenic heat flux can be much stronger (Sailor and Lu 2004). Ichinose et al (1999) showed that the anthropogenic heat flux in central Tokyo exceeded 400 W m −2 in the daytime, and the maximum value reached 1590 W m −2 in the early morning in winter.
Previous studies on the effects of anthropogenic heat flux on urban climate were typically conducted using weather and climate models. Salamanca et al (2014) quantified the impacts of anthropogenic heat flux by turning on/off air conditioning systems in the Weather Research and Forecasting (WRF) model coupled to a building energy model (BEM) and a multilayer building effect parameterization. Their results revealed that the heat emitted from air conditioning systems resulted in a 1 • C-1.5 • C temperature rise during summer nights over Phoenix, USA. Fan and Sailor (2005) incorporated an anthropogenic heating source term in the near-surface energy balance within the NCAR/PennState Fifth Generation Model. They found that the influence of anthropogenic heat flux on the urban climate of Philadelphia, USA was significant, particularly during nighttime and in winter, with near-surface air warming as large as 2 • C-3 • C. Similar results were also found in China (Feng et al 2012(Feng et al , 2014 and Australia (Ma et al 2017), where the temperature rise was more pronounced in winter than summer. In another numerical study conducted in a Japanese megacity (Keihanshin district), the results indicated that although the daytime anthropogenic heat flux was larger than the nighttime counterpart, the induced temperature rise was nearly threefold larger at night (Narumi et al 2009). Studies also revealed that the anthropogenic heating effects depended not only on the quantity of anthropogenic heat flux but also atmospheric stratification and orographic factors (Block et al 2004, Narumi et al 2009, Zhang et al 2016. Since it is obvious that the amount of warming induced by anthropogenic heating depends on the magnitude of the anthropogenic heat flux, it is perhaps more important to examine the ratio of the temperature increase to the amount of anthropogenic heat flux (∆T a /∆Q AH , where ∆ represents a change), much like the concept of climate sensitivity but at a local (urban) scale. In this sense, we treat the change in anthropogenic heat flux (∆Q AH ) as the climate forcing and the change in urban temperature (∆T a ) as the climate response. Table 1 provides a selected list of existing studies on the warming effect of anthropogenic heat flux. By normalizing the temperature increase by the magnitude of anthropogenic heat flux, a better consistency among different studies emerges, with the magnitude of ∆T a /∆Q AH being of the order of 0.01 K ( W m −2 ) −1 . This value is consistent with the findings of Kikegawa et al (2014), who carried out field campaigns based on meteorological measurements and monitoring of electricity demand, as well as numerical simulations with WRF (coupled with a multilayer urban canopy model and a BEM) in two major Japanese cities, Tokyo and Osaka, in July to August 2007. Their work suggested an afternoon sensitivity of 0.01 K ( W m −2 ) −1 based on observations and showed that the simulated results had the same order of magnitude. However, it is noteworthy to point out that the magnitude of ∆T a /∆Q AH from different studies (table 1) still varies by nearly two orders of magnitude [from 0.001 to 0.05 K ( W m −2 ) −1 ]. More importantly, the physical processes responsible for such variability remain elusive. Quantification of the key factors controlling the variability of ∆T a /∆Q AH frames the scope of this study.
To achieve this, we develop a forcing-feedback framework based on the energy budget of air within the urban canopy layer (UCL, i.e. the layer below the height of the main urban elements) and apply it to diagnosing ∆T a /∆Q AH simulated by the Community Land Model Urban (CLMU) over the CONUS. This region, characterized by a growing urban population and significant energy consumption, has not yet been thoroughly investigated in terms of the impact of anthropogenic heat flux. This study is organized as follows: section 2 describes the forcing-feedback framework and model experiments. Section 3 evaluates ∆T a /∆Q AH at the seasonal and diurnal scales. The key feedback mechanisms and the factors controlling the variability of ∆T a /∆Q AH are discussed in detail in this section. Finally, discussions and conclusions are presented in sections 4 and 5, respectively.

A forcing-feedback framework
We propose a forcing-feedback framework to diagnose the sensitivity of air temperature within the UCL, also called urban canopy air temperature hereafter, to anthropogenic heat flux based on the energy budget of air within the UCL (figure 1). This conceptualization of the UCL is consistent with the theoretical underpinning of nearly all single-layer urban canopy models (UCMs) in weather and climate modeling, including the CLMU to be used in this study (more details on the CLMU are presented later). Our starting point is that the UCL is our control volume (or system of interest) and is the direct recipient of anthropogenic heat flux (i.e. the forcing). At steady state, the energy budget of the air within the UCL can be written as where Q AH is the anthropogenic heat flux and R is the sum of heat fluxes other than the anthropogenic heat flux (we shall say more about R later). When the anthropogenic heat flux is altered by a certain amount (indicated by ∆), the energy balance of air within the UCL reaches a new equilibrium state where ∆Q AH can be interpreted as the added anthropogenic heat flux compared with the scenario without anthropogenic heat flux and ∆R is the total change of other heat fluxes in response to ∆Q AH . Changes in other heat fluxes (∆R) are often related to changes in the urban canopy air temperature (∆T a ). Denoting ∆R = λ all ∆T a , we can write the sensitivity of urban canopy air temperature to anthropogenic heat flux as where λ all is the sensitivity parameter (called the total sensitivity parameter in order to distinguish it from other feedback parameters introduced later). The sensitivity ∆T a /Q AH , which indicates how easily the urban canopy air temperature can be altered by a perturbation of anthropogenic heat flux, is thus equivalent to the negative reciprocal of the total sensitivity parameter (λ all ). If the absolute value of λ all is larger, the urban canopy air warming per unit increase of anthropogenic heat flux is weaker. Therefore, to understand ∆T a /Q AH , we need to examine λ all in the relation ∆R = λ all ∆T a . The air within the UCL receives convective heat fluxes from various urban surfaces and the overlying atmosphere. The sum of these heat fluxes (R) received by the air within the UCL can thus be written as where n is the number of urban surfaces (e.g. there are five urban surfaces in the CLMU that interact with the urban canopy air, including roof, previous ground, imperious ground, sun wall and shade wall), i refers to the ith urban surface, w i is the weight of the ith surface based on the corresponding area fraction (converted to per unit area of urban canyon floor in the horizontal direction), ρ is the air density (kg m −3 ), C p is the specific heat of air at constant pressure, assumed to be of a constant value of 1004.64 J kg −1 K −1 , c s is the heat conductance between the air within the UCL and the urban surface (called the surface-canopy air heat conductance, m s −1 ), c a is the heat conductance between the air within the UCL and the overlying atmosphere (called the atmosphere-canopy air heat conductance, m s −1 ), T s is the urban surface temperature (K) and θ atm is the atmospheric potential temperature (K). Here we have assumed that the heat conductances between the air within the UCL and different urban surfaces are identical, which is a common assumption made in the CLMU and many other single-layer UCMs. But this assumption can be relaxed by allowing c s to vary for different urban surfaces in future work. With equation (4), λ all can be written as the sum of the direct effect and feedbacks. Using the chain rule on equation (4) yields where the partial and total derivatives are denoted by ∂ and ∆, respectively. In this equation, λ 0 is the baseline sensitivity parameter, representing the direct effect of anthropogenic heat flux on urban canopy air temperature, with everything else (e.g. surface temperature, atmosphere-canopy air heat conductance, etc) held the same. Other λ parameters represent different feedback processes: λ 1 refers to the strength of feedback from changes in surface temperatures; λ 2 is the feedback parameter for changes in atmosphere-canopy air heat conductance; λ 3 is the feedback parameter for changes in surface-canopy air heat conductance; λ 4 is the parameter for atmospheric feedback. A positive (or negative) feedback means that the process leads to an amplification (or dampening) of the direct effect of anthropogenic heat flux on urban canopy air temperature.
Combining equations (4) and (5), the baseline sensitivity parameter and feedback parameters can be derived as Equations (1)-(10) constitute our forcingfeedback framework for diagnosing the sensitivity of urban canopy air temperature to anthropogenic heat flux. The aim of the proposed forcing-feedback framework is not to predict ∆T a /∆Q AH but to provide a diagnostic tool for quantifying the strengths of direct effects and feedback processes. In this study, the inputs for this framework are the simulated results from the CLMU. However, this framework is not limited to the CLMU and can be applied to diagnosing outputs from other UCMs.

The CLMU model and the numerical experiment design
The CLMU model represents urban parameterization within the Community Land Model (CLM), which is the land component of the Community Earth System Model (CESM) (Danabasoglu et al 2020). In this study, the most recent released version of the CLM (CLM5) within the framework of CESM version 2 (CESM2) is used. Within each land grid cell, CLM5 can have multiple land units including vegetated, crop, urban, glacier and lakes. For each urban land unit, three urban categories (tall building district, high density and medium density) are allowed.
In CLMU, the urban canyon system consists of five surfaces: roofs, sunlit and shaded walls, impervious and pervious floors. The energy and water fluxes from each urban surface interact with the canopy air (see figure 1). A more detailed description of CLMU including the main urban parameters can be found elsewhere Feddema 2020). The CLMU input data are supplied by a global dataset (Jackson et al 2010). The model has been widely used to study urban energy and water fluxes, as well as surface and air temperatures (Oleson et al 2008a, 2008b, Demuzere et al 2013, Karsisto et al 2016, Oleson and Feddema 2020. In this study, we use an improved CLMU that includes parameterizations of urban heat mitigation strategies (e.g. cool roofs and green roofs) that have been proposed and validated in our previous work (Wang et al , 2021, although these new features are not used in this study. We run CLM5 in an offline mode (i.e. forced by meteorological data) at a 1/8 degree spatial resolution over the CONUS and at an hourly time step (Wang and Li 2022). The hourly meteorological forcing data are from the North America Land Data Assimilation System phase II (NLDAS2) dataset (Xia et al 2012). The model is first spun up for 84 years by recycling the 1979-1999 NLDAS2 forcing four times. Four sets of numerical experiments are then conducted from 1979 to 1999 using the same initial condition obtained from the spin-up run (table S1). These four numerical experiments are designed to quantify how the urban canopy air temperature (figure 1) responds to a prescribed increase anthropogenic heat flux. In the control (CTL) experiment, no anthropogenic heat flux is added to the urban canopy air heat budget, and the simulated canopy air temperature is denoted as T a,0 . In the first sensitivity experiment (AH1), we add 1 W m −2 of anthropogenic heat flux into the urban canopy air heat budget at each time step and compute a new canopy air temperature (hereafter T a,1 ). Therefore, the difference between T a,1 and T a,0 (which is numerically equivalent to ∆T a /∆Q AH given that the added anthropogenic heat flux is 1 W m −2 ) is the total impact of 1 W m −2 of anthropogenic heat flux, which includes both the direct effect and the feedbacks. In another two sensitivity experiments AH10 and AH100, the added anthropogenic heat flux is 10 W m −2 and 100 W m −2 , respectively. We denote the simulated urban canopy air temperatures in these two experiments as T a,10 and T a,100 , respectively. The sensitivity ∆T a /∆Q AH is thus calculated as (T a,10 − T a,0 ) /10 and (T a,100 − T a,0 ) /100, respectively (see table S1). These two sensitivity experiments (AH10 and AH100) are designed to quantify whether the sensitivity ∆T a /∆Q AH is influenced by the magnitude of anthropogenic heat flux due to nonlinearity in the feedback processes. We choose these values (1, 10, 100 W m −2 ) to cover a wide but reasonable range of anthropogenic heat fluxes.
We should emphasize that the added anthropogenic heat flux in all our experiments is prescribed, not computed by the BEM in CLMU (Oleson et al 2011, Demuzere et al 2013. We prescribe the added anthropogenic heat flux because we are mostly interested in the sensitivity of urban canopy air temperature to anthropogenic heat flux, not what processes generate the anthropogenic heat flux. Moreover, when the BEM in CLMU is used, the generated anthropogenic heat flux is added to the pervious and impervious surface energy budgets, which seems unphysical and is avoided in our study. Another way of interpreting our results is that they represent the sensitivity of urban canopy air temperature to anthropogenic heat fluxes from non-building (e.g. transportation) sectors with magnitudes of 1, 10 and 100 W m −2 .
For all simulations, we output the hourly urban canopy air temperature, the temperatures of different urban surfaces (i.e. roofs, walls and canyon floors), the atmospheric potential temperature, the surfacecanopy air heat conductance (c s ), as well as the atmosphere-canopy air heat conductance (c a ). Note that the outputted c s and c a are computed internally via their parameterizations in CLMU. These hourly outputs are then used in the forcing-feedback framework described in section 2.1. Specifically, we compute the hourly sensitivity parameters based on equations (6)-(9). Given that we do not have atmospheric feedbacks in our simulations, λ 4 = 0. With the sensitivity parameters calculated using equations (6)-(9), the total sensitivity ∆T a /∆Q AH can be diagnosed using equations (3) and (5). The diagnosed ∆T a /∆Q AH is then compared with the directly computed ∆T a /∆Q AH mentioned above (e.g. for AH1 the directly computed ∆T a /∆Q AH is simply T a,1 − T a,0 ). We average the hourly results over 20 years from 1980 to 1999.
Before we move to the results section, it is informative to briefly discuss the physics behind c a and c s and their parameterizations in CLMU, as they are key parameters in the forcing-feedback framework (see, e.g., equation (6)). Physically, c a (c s ) represents the efficiency of convective heat transfer between the overlying atmosphere (the urban surfaces) and the canopy air. Given that the flow within the UCL is turbulent, both c a and c s are strongly affected by shear and buoyancy, the two main sources of turbulence kinetic energy. However, c a and c s are fundamentally different because they represent the convective heat transfer efficiencies across different levels. In terms of their parameterizations in CLMU, c a is parameterized through the classic Monin-Obukhov similarity theory (Oleson et al 2008a). Hence, c a is strongly affected by atmospheric stratification. However, c s is parameterized as a function of wind speed alone in the urban canyon (Oleson et al 2008a) and is thus much less affected by atmospheric stratification than c a .

Sensitivity of urban canopy air temperature to anthropogenic heat flux (∆T a /∆Q AH ) and the associated feedback parameters
We first present the sensitivity ∆T a /∆Q AH simulated by CLMU in summer (June-August, or JJA) and winter (December-February, or DJF) seasons (figure 2). The results shown here have been averaged over 20 years (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999) and are based on the AH1 experiment where the added anthropogenic heat flux is 1 W m −2 . The effect of increasing the magnitude of anthropogenic heat flux will be discussed in section 3.4. In summer, the median value of ∆T a /∆Q AH is around 0.01 K ( W m −2 ) −1 , broadly comparable with previous studies presented in table 1. Here the median values are shown to minimize the influence of outliers. In winter, ∆T a /∆Q AH becomes stronger, with the median value increased by about 20%. In some cities in the southwestern USA (e.g. Los Angeles and Phoenix), the winter values of ∆T a /∆Q AH even reach 0.03 K ( W m −2 ) −1 . To understand the directly computed ∆T a /∆Q AH from CLMU simulation results, we employ the forcing-feedback framework described in section 2.1. The total sensitivity diagnosed from this framework (i.e. using equations (5)-(9)) matches very well with the directly computed ∆T a /∆Q AH (figure 2), with spatial correlation coefficients larger than 0.99. These results give us confidence to use the forcing-feedback framework to interpret ∆T a /∆Q AH .
Based on the forcing-feedback framework, the total sensitivity parameter (λ all ) can be decomposed into the sum of the baseline sensitivity parameter (λ 0 ) and the feedback parameters (λ 1 , λ 2 and λ 3 ). We find that the magnitude of λ all is almost identical to the magnitude of λ 0 (the spatial median value is −122 W m −2 K −1 for both λ all and λ 0 ) in summer (figure 3). This is because the sum of the three feedback parameters (λ 1 + λ 2 + λ 3 ) is very small, with the positive feedbacks and negative feedbacks nearly canceling each other. The positive feedback is mainly from changes in surface temperature (λ 1 , with a median value of 24 m −2 K −1 ). This is expected as increases in surface temperature due to the added anthropogenic heat flux can in turn amplify the urban canopy air warming. On the other hand, the negative feedback is mainly the result of changes in atmosphere-canopy air heat conductance (λ 2 , with a median value of −25 m −2 K −1 ). As the anthropogenic heat flux is added, the atmospheric stratification is altered (i.e. relatively more unstable), resulting in increased atmosphere-canopy air heat conductance (c a ). This in turn leads to an increase in heat transfer into the overlying atmosphere and a dampening of the urban canopy air warming signal. The feedback from changes in surface-canopy air heat conductance (λ 3 , with a median value of 1 m −2 K −1 ) is much weaker than the other two feedback processes. This can be explained by the parameterization of surfacecanopy air heat conductance (c s ) in CLMU, which is only dependent on the wind speed in the UCL and is thus a much weaker function of atmospheric stratification than the atmosphere-canopy air heat conductance (c a ). In winter (figure 3), the negative feedback from atmosphere-canopy air heat conductance (λ 2 ) decreases in magnitude by 11 W m −2 K −1 (in terms of median value) when compared with its summer counterpart (see also figure S1 for a comparison between summer and winter results). Namely, λ 2 becomes less negative, implying that the negative feedback from atmosphere-canopy air heat conductance (c a ) is weakened. Unlike the reduced magnitude of λ 2 in winter, the winter-summer differences in λ 1 and λ 3 are much smaller (about 1 W m −2 K −1 in terms of median values) and almost negligible. As a result, the sum of feedbacks (λ 1 + λ 2 + λ 3 ) becomes positive in winter (compared with nearly zero in summer). The absolute value of the total sensitivity parameter (λ all ) therefore decreases, which further leads to an increase in ∆T a /∆Q AH . The weakened negative feedback from atmosphere-canopy air heat conductance (c a ) explains why stronger canopy air warming is observed in winter than in summer with the same amount of anthropogenic heat flux (figure 2), a typical result in the literature. Figure 2 shows the strong spatial variabilities in ∆T a /∆Q AH . To understand these spatial variabilities, we first note that the spatial pattern of the baseline sensitivity parameter λ 0 is very close to that of λ all , with spatial correlation coefficients of 0.77 and 0.95 in summer and winter, respectively. Therefore, the spatial variability of λ 0 largely determines the spatial variability of ∆T a /∆Q AH . From equation (6), λ 0 is proportional to the sum of atmosphere-canopy air heat conductance (c a ) and surface-canopy air heat conductance (c s ). We find that c s is less than 20% of c a and shows little spatial variability (not shown). As a result, one would expect that the spatial variability of λ 0 is mainly controlled by the spatial variability of c a . This is indeed the case. We find that the spatial correlation coefficients between λ 0 and c a are very strong (−0.87 and −0.98 in summer and winter, Figure 3. The sensitivity and feedback parameters: (a), (b) the total sensitivity parameter (λ all ); (c), (d) the baseline sensitivity parameter (λ0); (e)-(g) the feedback parameter for surface temperature (λ1) ((e), (f)), heat conductance between the canopy air and overlying atmosphere (λ2) ((g), (h)), and heat conductance between the canopy air and urban surfaces (λ3) ((i), (j)). Parts respectively). The negative correlations are understandable since physically the atmosphere-canopy air heat conductance (c a ) indicates how strongly the air within the UCL communicates with the overlying atmosphere in terms of convective heat transfer. In places with larger (smaller) c a , it is easier (more difficult) to transfer heat from the UCL to the overlying atmosphere, and thus the canopy air warming signal is weaker (stronger) with the same amount of anthropogenic heat flux.

Diurnal variation of ∆T a /∆Q AH and its controlling factors
We further analyze the diurnal variation of ∆T a /∆Q AH . To do this, we select four metropolitan cities that have widely different climates and geographical locations (San Francisco, Boston, Chicago and Houston), instead of presenting averaged results over the CONUS.
In summer ( figure 4(a)), all four cities experience a higher ∆T a /∆Q AH in the early morning than at other times. The morning peak of ∆T a /∆AH is around 0.038 K ( W m −2 ) −1 in Houston, followed by San Francisco, Boston and Chicago. In the afternoon, the sensitivity in all cities is close to 0.01 K ( W m −2 ) −1 , which is consistent with the findings of Kikegawa et al (2014) that also suggested a summer afternoon sensitivity of 0.01 K ( W m −2 ) −1 . In contrast, there exist large differences in the diurnal variation of ∆T a /∆Q AH in winter ( figure 4(b)). The sensitivity ∆T a /∆Q AH in Boston and Chicago is around 0.01 K ( W m −2 ) −1 throughout the day with small diurnal variations, while the diurnal variations of ∆T a /∆Q AH in Houston and San Francisco are strong, with much larger nighttime values than daytime values. San Francisco has the largest sensitivity in winter among the four cities, with a peak value of 0.036 K ( W m −2 ) −1 . According to the forcing-feedback framework, the diurnal variations of ∆T a /∆Q AH are linked to the diurnal variations of feedback parameters, including the baseline sensitivity parameter (λ 0 ). As shown in figures 4(c)-(l), λ 0 and, to a lesser extent, λ 2 exhibit diurnal variations that resemble those of λ all , implying that the diurnal variations of ∆T a /∆Q AH are controlled by processes encoded in λ 0 (equation (6)) and, to a lesser extent, λ 2 (equation (8)). Close inspection of equations (6) and (8) indicates that a common process in equations (6) and (8) is the atmospherecanopy air heat conductance (c a ), suggesting that the diurnal variation of c a (and ∆c a ) are the key to understanding the diurnal variation of ∆T a /∆Q AH .
The c a is controlled by shear-and buoyancygenerated turbulence and thus is strongly affected by atmospheric stratification. In winter, the air within the UCL experiences more stable conditions at night, and hence c a is smaller, λ 0 is less negative (figure 4(f)) and ∆T a /∆Q AH is larger ( figure 4(b)) than their daytime counterparts, assuming that the shear is the same between daytime and nighttime. In summer, the accumulation of stable stratification throughout the night reduces c a (leading to a less negative λ 0 , figure 4(e)) and increases ∆T a /∆Q AH (figure 4(a)). After sunrise, the stratification transitions from stable to unstable, which increases c a , causes a more negative λ 0 and reduces ∆T a /∆Q AH . These two processes yield a morning peak of ∆T a /∆Q AH , as observed in figure 4(a). Shear also plays an important role. For example, the stronger winds in Boston and Chicago in winter are likely to cause larger shear, leading to a larger c a and smaller ∆T a /∆Q AH , when compared with Houston and San Francisco ( figure 4(b)).

Nonlinear response of ∆T a to ∆Q AH
The above results are from the AH1 experiment, which adds 1 W m −2 of anthropogenic heat flux into the UCL. We also conduct experiments to investigate how the urban canopy air temperature responds to different amounts of anthropogenic heat flux. The aim of these experiments is to test whether any of the feedbacks scale nonlinearly with ∆Q AH , thereby creating nonlinear responses of ∆T a to ∆Q AH . Note that the baseline sensitivity parameter (λ 0 , see equation (6)) does not change with the magnitude of anthropogenic heat flux. Thus, any nonlinear response must stem from the feedback processes. Figure 5 presents the relative changes in ∆T a /∆Q AH and feedback parameters (λ 1 , λ 2 , λ 3 ) by comparing AH10 and AH100 with AH1 (i.e. the results of AH10 and AH100 minus the results of AH1 and then normalized by the results of AH1). The relative changes in ∆T a /∆Q AH are all negative, implying that the sensitivity becomes smaller as the magnitude of anthropogenic heat flux increases. The relative changes between AH100 and AH1 in terms of ∆T a /∆Q AH have median values of −27% and −35% in summer and winter, respectively. This suggests that ∆T a does respond nonlinearly to ∆Q AH . Here we should stress that this result does not mean that changes in urban canopy air temperature ∆T a become smaller as the magnitude of anthropogenic heat flux increases. It is rather ∆T a /∆Q AH that reduces as the magnitude of the anthropogenic heat flux increases.
The relative changes in feedback parameters suggest that the nonlinear response of urban canopy air warming to the addition of anthropogenic heat flux is mostly due to decreases in λ 2 (i.e. λ 2 becomes more negative) as ∆Q AH increases (figure 5). For example, the differences between AH100 and AH1 in terms of λ 2 give median values of −13% and −28% in summer and winter, respectively. As alluded to earlier in section 3.1, λ 2 is associated with changes in the atmosphere-canopy air heat conductance (∆c a ). These results imply that with a larger ∆Q AH , the increase in c a is stronger, leading to a more negative λ 2 and a weaker ∆T a /∆Q AH . Therefore, the nonlinear response of ∆T a to ∆Q AH is traced to the role of c a .

Discussion
This study has several implications that are important to appreciate. First, we emphasize that it is equally important to study the sensitivity (∆T a /∆Q AH ) in addition to the forcing magnitude (∆Q AH ). The sensitivity is the ratio of the response (∆T a ) to the forcing and is a much better constrained quantity than the response itself, as can be seen from table 1. Second, the forcing-feedback framework further allows us to understand why many previous studies reported a ∆T a /∆Q AH value of approximately 0.01 K ( W m −2 ) −1 . Without considering any feedbacks and any role for c s (both are reasonably good assumptions), the baseline sensitivity Third, the forcing-feedback framework allows us to quantify the contributions of various physical processes to the spatiotemporal variability of ∆T a /∆Q AH . Our results demonstrate that the atmosphere-canopy air heat conductance (c a ) plays a key role in controlling the spatiotemporal variations of ∆T a /∆Q AH , as well as the nonlinear response of ∆T a to ∆Q AH . Hence, it is critical for UCMs to accurately represent the convective heat transfer between the canopy air and the overlying atmosphere, among other things. Currently, Monin-Obukhov similarity theory remains the workhorse model for parameterizing c a in UCMs due to its popularity and parsimony (e.g. in the CLMU, see Oleson et al 2008a), even though urban areas are not homogeneous and thus Monin-Obukhov similarity theory does not strictly apply (Garratt 1994). It remains unclear whether Monin-Obukhov similarity theory combined with urban roughness lengths is sufficient for parameterizing c a over urban areas or if new theories accounting for the effects of urban canopies (e.g. similar to the work by Finnigan (2007, 2008), see also Bonan et al (2018)) are needed. Furthermore, in this context nearly all UCMs assume that turbulent transport is the only process that needs to be parameterized. However, dispersive transport might become relevant over areas with large variations of building heights (Akinlabi et al 2022). Addressing these questions is outside the scope of this study but is strongly needed.
There are also limitations of this work that need to be pointed out. First, we only evaluate the feedback processes within the UCL. Quantifying the role of atmospheric feedback (λ 4 ) and how it is scale-dependent (Li and Wang 2019) is left for future work. Second, while we highlight the key role played by the atmosphere-canopy air heat conductance (c a ), diagnosing the physical processes as well as urban morphological parameters that give rise to the spatiotemporal variability of c a (e.g. diagnosing the differences between different cities in figure 4) remains to be conducted. Within the confines of Monin-Obukhov similarity theory, c a is affected by sheargenerated and buoyancy-generated turbulence and is a function of mean wind speed, roughness lengths (both momentum and thermal roughness lengths) and stability parameters. The momentum roughness length is further a complex function of building height and canyon geometry. Understanding the spatiotemporal variability of c a and its relation to these underlying factors is beyond the scope of this study. Third, this study does not prescribe spatially and temporally varying anthropogenic heat flux. This is justified by the focus of this work on the sensitivity (∆T a /∆Q AH ) instead of the response (∆T a ). ∆T a can be viewed as the product of ∆T a /∆Q AH and the forcing (∆Q AH ). Thus, the spatiotemporal variability of the temperature response is further complicated by the spatiotemporal variability of the forcing. Studies aiming to quantify the temperature response should also address the variability of the forcing.

Conclusion
Anthropogenic heat flux is an important controlling factor of the urban thermal environment. Although many studies have investigated the impacts of anthropogenic heat flux, the key factors controlling the magnitude of the sensitivity of urban air temperature to anthropogenic heat flux (∆T a /∆Q AH ) and its spatial and temporal patterns remain elusive. In this study, we develop a forcing-feedback framework based on the energy balance of air within the UCL and apply the framework to diagnosing simulated ∆T a /∆Q AH over the CONUS by a numerical model. Within the forcing-feedback framework, ∆T a /∆Q AH is decomposed into the direct effect of Q AH on T a , as well as feedbacks through changes in the surface temperature (T s ), the atmosphere-canopy air heat conductance (c a ), and the surface-canopy air heat conductance (c s ). This forcing-feedback framework allows us, for the first time, to understand the contributions of physical processes within the UCL to ∆T a /∆Q AH and the spatiotemporal variability of ∆T a /∆Q AH in a quantitative manner.
Our study first examines the seasonal variation of ∆T a /∆Q AH . In summer, the positive feedback (mainly from changes in surface temperature, represented by λ 1 ) is nearly canceled by the negative feedback (mainly from changes in atmosphere-canopy air heat conductance c a , represented by λ 2 ). As a result, ∆T a /∆Q AH is dominated by the direct effect (represented by λ 0 ). In winter, the negative feedback from c a (represented by λ 2 ) weakens, leading to a stronger ∆T a /∆Q AH . We also investigate the diurnal variations of ∆T a /∆Q AH . The results show that the diurnal variations of ∆T a /∆Q AH are mostly controlled by the diurnal variations in λ 0 , and to a lesser extent, λ 2 , both of which are strongly related to the diurnal variations of c a (and ∆c a ). Hence, it can be summarized that the temporal (both seasonal and diurnal) dynamics of ∆T a /∆Q AH are mostly controlled by those of c a . We also find that the spatial variability of ∆T a /∆Q AH over the CONUS is mainly determined by the direct effect (λ 0 ). Since λ 0 is proportional to the sum of c a and c s , and c s shows little spatial variability, the spatial variability of ∆T a /∆Q AH is dominated by the spatial variability of c a . We further examine the nonlinearity in the response of ∆T a to ∆Q AH by varying the magnitude of ∆Q AH . The nonlinear response of ∆T a to ∆Q AH stems mostly from the feedback process associated with changes in atmosphere-canopy air heat conductance (c a ). Our framework provides a tool for studying the feedback mechanisms that are important for understanding the sensitivity of urban canopy air temperature to anthropogenic heat flux.