The Nonlinear Radiative Feedback Effects in the Arctic Warming

The analysis of radiative feedbacks requires the separation and quantification of the radiative contributions of different feedback variables, such as atmospheric temperature, water vapor, surface albedo, cloud, etc. It has been a challenge to include the nonlinear radiative effects of these variables in the feedback analysis. For instance, the kernel method that is widely used in the literature assumes linearity and completely neglects the nonlinear effects. Nonlinear effects may arise from the nonlinear dependency of radiation on each of the feedback variables, especially when the change in them is of large magnitude such as in the case of the Arctic climate change. Nonlinear effects may also arise from the coupling between different feedback variables, which often occurs as feedback variables including temperature, humidity and cloud tend to vary in a coherent manner. In this paper, we use brute-force radiation model calculations to quantify both univariate and multivariate nonlinear feedback effects and provide a qualitative explanation of their causes based on simple analytical models. We identify these prominent nonlinear effects in the CO2-driven Arctic climate change: 1) the univariate nonlinear effect in the surface albedo feedback, which results from a nonlinear dependency of planetary albedo on the surface albedo, which causes the linear kernel method to overestimate the univariate surface albedo feedback; 2) the coupling effect between surface albedo and cloud, which offsets the univariate surface albedo feedback; 3) the coupling effect between atmospheric temperature and cloud, which offsets the very strong univariate temperature feedback. These results illustrate the hidden biases in the linear feedback analysis methods and highlight the need for nonlinear methods in feedback quantification.


INTRODUCTION
Radiative forcing and feedbacks strongly influence the Arctic climate. The warming in the Arctic has occurred in a faster pace than the global average, due to greenhouse gas forcing and amplifying feedbacks (Stocker et al., 2013). It requires accurate quantification of the radiative effects of associated feedback variables (surface albedo, atmospheric temperature, water vapor, cloud, etc.,) in order to ascertain their contributions to the climate change of interest. For instance, based on the energy budget balance with regard to the Top-of-Atmosphere (TOA), surface or atmospheric budget and assuming the warming induced thermal radiation (Planck effect) balances the radiation changes caused by feedbacks, one can infer how much global or regional warming, e.g., the Arctic warming amplification, can be attributed to individual feedbacks (Held and Soden 2000;Lu and Cai 2009;Pithan and Mauritsen 2014).
Often assumed in feedback analysis is linear additivity of the radiative effects of different feedback variables. For instance, the widely adopted kernel method (Soden and Held 2006) measures the radiation change caused by a feedback variable (X) by multiplying a pre-calculated radiative kernel ( zR zX ) with the climate response (dX). Due to its simple concept and computational efficiency, a large number of studies have been conducted using this method (e.g., Soden and Held 2006;Zelinka et al., 2012;Vial et al., 2013;Zhang and Huang, 2014) and precomputed kernels based on different atmospheric datasets, including climate models, reanalyses and satellite data (e.g., Soden et al., 2008;Yue et al., 2016;Huang et al., 2017).
The nonlinear effects, however, are often too large to ignore. When individual feedback terms are independently measured, such as the non-cloud feedbacks in the clear-sky case in the kernel method, the ignored nonlinear effects may lead to a non-closure of the radiation budget, i.e., the sum of the individual terms cannot reproduce the overall radiation change (e.g., Huang 2013;Vial et al., 2013). In the Arctic, where climate perturbations are of large magnitudes, e.g., in the case of sea ice melt, the non-closure issue is especially noticeable (e.g., Shell et al., 2008;Block and Mauritsen 2013;Zhu et al., 2019). Besides the large perturbations in surface albedo, the Arctic is also noted for its strong and unique lapse rate (e.g., Pithan and Mauritsen 2014) and cloud (e.g., Kato et al., 2006) feedbacks. It should be noted that although some methods exhibit a seemingly good radiation closure, the nonlinear effects are not treated but hidden in the feedback term(s) measured as a residual, e.g., the cloud feedback term in the typical kernel method, including the adjusted cloud radiative forcing (aCRF) technique Soden et al., 2008).
Although the existence of the nonlinear effects has been recognized (e.g., Zhang et al., 1994;, their impacts were seldom isolated and quantified. Some recent works have specifically addressed the nonlinearity issue in the radiative feedback analysis. Zhu et al. (2019) for the first time used a neural network model (a nonlinear diagnostic method without linearity assumption) to assess the radiative feedbacks and identified a few strong nonlinear effects, including a strong cloud-water vapor coupling effect in the tropical climate variations and a strong nonlinear dependence of radiation flux on the surface albedo. Using Partial-Radiative-Perturbation (PRP) experiments and brute-force radiation model-based computations, Huang and Huang (2021) verified the cloud-water vapor coupling effect and offered an analytic estimation of this effect on the longwave radiation. Shakirova and Huang (2021) advanced the neural network model of Zhu et al. and demonstrated its advantages particularly for quantifying the albedo feedback.
In this paper, we aim to give an overview of the nonlinear radiative feedback effects in Arctic climate change. Based on a heuristic climate change scenario of broad interest: the abrupt quadrupling of atmospheric CO 2 (4xCO 2 ), we investigate how the nonlinear radiative effects arise from the univariate and multivariate variations of the feedback variables, such as atmospheric and surface temperature (t), water vapor (q), surface albedo (a) and cloud (c), and measure how the nonlinear effects compare to the linear effects in terms of magnitude and pattern. We note that in this paper we are not concerned with how the changes in these variables are resulted, which if nonlinearly related to the surface warming may also cause nonlinearity in climate feedbacks, but focus on how their changes, as projected by the GCM, lead to nonlinear changes in the TOA longwave (LW) and shortwave (SW) radiation energy fluxes. In the following sections, we will define, demonstrate and discuss the various feedback effects of interest in order.

METHOD: FEEDBACK DEFINITIONS
Here, we define a radiative feedback as the (partial) radiation change, in the units of W m −2 , due to one or multiple feedback variables. This should be distinguished from a feedback parameter, which is normalized by surface temperature change and is in the units of W m −2 K −1 .
Consider the radiation field of interest, e.g., the TOA or surface radiation flux, as a function of the feedback variables: R R(x, y, z), where the letters (x, y, z) are generic notations of the feedback variables. The total radiation change in a given climate change scenario can thus be expressed by a Taylor series as where the subscripts 1 and 2 denote two different climate states, e.g., those before and after quadrupling CO 2 (noted as 1xCO 2 and 4xCO 2 , respectively from now on); such terms as Δx x 2 − x 1 denote the climate responses. The terms on the righthand side of Eq. (1) illustrate three types of radiative effects that we aim to elucidate here: 1) The univariate linear effects, such as zR zx Δx, which we denote as ΔR x ; 2) The univariate nonlinear effects, such as 1 2 z 2 R zx 2 (Δx) 2 , which we denote as ΔR xx ; 3) The multivariate nonlinear effects, such as z 2 R zxzy ΔxΔy, which we denote as ΔR xy .

zR zx
Δx + 1 2 August 2021 | Volume 9 | Article 693779 Equation (2) is written following the PRP concept (Wetherald and Manabe 1988) and measures the univariate feedbacks by simply evaluating the radiation flux twice: first with an unperturbed profile, R(x 1 , y 1 , ...), and then perturbing x only, R(x 2 , y 1 , ...). Note that other unperturbed independent variables than y are omitted in these expressions. If not otherwise stated, unspecified independent variables all take the unperturbed values when the radiation fluxes are evaluated in the following. The evaluation of the radiation fluxes can be done using a physical model, i.e., a radiative transfer model (RTM) , or a statistical model, e.g., a neural network model that emulate the radiation fluxes (Zhu et al., 2019). Similarly, a bivariate feedback can be expressed as: From Eq. 2 and Eq. 3, the bivariate coupling effect ΔR xy can be obtained as Based on the above equations and following Huang and Huang (2021), we evaluate the radiation fluxes and isolate the respective feedback effects (ΔR x , ΔR xx , ΔR xy , etc.), using the Rapid Radiative Transfer Model (RRTM, Mlawer et al., 1997). Using this RTM, the TOA and surface radiation fluxes are computed offline from instantaneous atmospheric profiles generated by the Community Earth System Model, CESM1.2, in a quadrupling CO 2 experiment (Wang and Huang, 2020). More details of the flux computation can be found in Huang and Huang (2021); we note that the radiative transfer computations, including the PRP computations, are based on instantaneous (3-hourly, as opposed to monthly mean) profiles at the original horizontal resolutions (1.9°× 2.5°) of the CESM and then averaged monthly or annually in all the results presented in the following section.
Note that we use a one-sided PRP, starting with the unperturbed climate and then prescribing the change(s) in the variables of interest, to define feedbacks, i.e., how much radiation change is caused by the change of the feedback variable(s) of concern. Some studies opt to use two-sided perturbations (e.g., . In contrast to Eq. (2), one may evaluate the feedback of x as which, as shown in the above expansion, effectively includes nonlinear coupling effects such as 1 2 z 2 R zxzy ΔxΔy in δR (x) . When individual feedbacks are evaluated this way, their sum can better reproduce the overall radiation change, i.e., achieving a better radiation closure. However, it should be noted these "individual" feedbacks δR (x) differ from the univariate feedbacks ΔR (x) as δR (x) contains coupling effects. To disclose these coupling effects, we adopt the one-sided formulations as exemplified by Eq. 2, Eq. 3, and Eq. 4

RESULTS
In this paper, we use the climate change in an abrupt 4xCO 2 experiment of CESM (Wang and Huang, 2020) to provide a context for examining the linear and nonlinear radiative feedbacks. As illustrated in Figure 1 for a few selected variables, this scenario represents strong perturbations in the Arctic climate, including reduction in surface albedo due to sea ice melt, surface warming and atmospheric moistening. The feedback quantifications presented in the following are based on the two months exemplified in Figure 1 if not otherwise noted.

Univariate Linear Effects
The univariate linear effect ΔR x can be measured, following its definition, by multiplying the radiative linear sensitivity kernel K x zR zx with the climate response Δx: R x K x Δx. This is the core idea of the kernel method (Soden et al., 2008). The kernels are usually pre-computed, again, following the PRP idea, by prescribing small perturbations to the individual variables, e.g., 1-K in atmospheric and surface temperatures, several percent change in water vapor concentration, or 0.01 increment of surface albedo (e.g., Shell et al., 2008): so that the kernel method is in essence to scale up the radiation change due to an infinitesimal (small) perturbation, R 0 , to estimate the radiation change due to a finite (large) perturbation: It is worth noting that when defining and applying the kernels, it is advisable to choose a scaling scheme appropriate to the radiation dependency on the feedback variable. For instance, in the case of such greenhouse gases as carbon dioxide and water vapor, their radiative effects are logarithmically dependent on their concentrations (e.g., Bani Shahabadi and Huang., 2014), so that it is common to define the water vapor kernel with respect to the change in the logarithm of the specific humidity, i.e., K q ΔR 0 Δ(ln(q)) 0 , or, in an approximate form, using the fractional change in q in the denominator, i.e., K q ΔR 0 (Δq/q) 0 . As shown by Figure 2, using the logarithmic scaling scheme, the radiation change caused by water vapor perturbations, even when the perturbations are of large magnitudes (about 200% increase of TCWV), can be well approximated according to Eq. (7). Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 693779 The validations against RTM-computed truth in Figures 2, 3 show that the non-cloud univariate radiative feedbacks ΔR (x) , even in the case of large climate perturbations, can be reasonably approximated by the linear term ΔR x . This is the basis of the kernel method (Soden et al., 2008). Nevertheless, the biases, i.e., the univariate nonlinear effects ΔR xx may amount to nonnegligible extent: in the case of the univariate water vapor feedback and the surface albedo feedback, the biases can amount to more than 10% in terms of Arctic mean (averaged over 70-90°N).
It should be cautioned that the kernel itself has a dependency on the atmospheric conditions (x, y, z, . . . ) and such dependency should be recognized when interpreting the kernel-diagnosed feedbacks. The kernels appropriate to evaluating the linear feedbacks in Eq. 1, Eq. 2, and Eq. 3, and used to compute ΔR q in Figure 2 and ΔR a in Figure 3, are computed from the GCM atmopsheric profiles in the unperturbed (1xCO2) climate, i.e., K x zR zx x 1 ,y 1 ,z 1 , where (x 1 , y 1 , z 1 ) denotes the 1xCO 2 climate. If computed from different atmospheric states, the kernel values may quantitatively differ (e.g., see the kernel comparisons in Huang et al., 2017;Smith et al., 2020 and others). If a kernel computed from a different atmospheric state (x 1 ′ , y 1 ′ , z 1 ′ ) is used to measure the linear feedback, a bias is resulted: Equation (8) shows that the bias can be considered one type of nonlinear effect in that it, like the nonlinear effects analyzed below, results from the nonlinear dependency of the radiation on the feedback variables (e.g., z 2 R zxzy ). For the simplicity of the expressions, we omit the notation (...)| x 1 ,y 1 ,z 1 in the following, where the conditioned states can be inferred from the context.
The magnitude of the feedback bias caused by kernel bias is proportional to the discrepancies in the atmospheric states (dx ′ , dy ′ , and dz ′ in Eq. (8)), which may introduce noticeable quantitative differences in the kernels. For example, one may see from Fig. S3 of Huang et al. (2017), as well as the discussions of Sanderson and Shell (2012), the temperature kernel discrepancies due to the ubiquitous discrepancies in the cloud distribution in different atmospheric datasets used for kernel computation. Here, as a sanity check, we recalculate  residual (A-B), i.e., the univariate nonlinear feedback, ΔR qq . Shown here is the clear-sky TOA LW radiation flux change due to the water vapor change in January in the 4xCO 2 experiment (as illustrated in Figure 1E). The kernel used in (B) is computed from the GCM instantaneous atmospheric profiles of the same month and is not subject to the bias discussed in Eq. 8.  Figure 1B).
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 693779 the Arctic mean clear-sky water vapor feedback shown in Figure 2 using the clear-sky kernels of Huang et al. (2017) and obtain an Arctic mean ΔR q of 5.79 W m −2 ( Figure 4). This, compared to the RTM-computed truth value of 3.05 W m −2 , representes a 90% bias and is much larger than the bias when the correct (GCM 1xCO 2 climate based) kernels are used (3.47 W m −2 , as shown by Figure 2B). The results indicates that, contradictory to common belief, there may be large biases, especially in regional (e.g., Arctic) feedbacks, resulting from kernel biases. Lastly, we note that cloud feedback is difficult, if not impossible, to be approximated by linear kernels. This is because cloud variations involve multiple radiative properties, including cloud fraction, droplet concentration and size distribution, etc., each of which may experience large, discrete perturbations and strongly affect the radiative sensitivity to each other. The cloud radiative effects measured in the cloud property histogram method (Zelinka et al., 2012) illustrate how the radiative sensitivity to cloud varies strongly with the cloud properties. Among other issues, a notable challenge is the vertical masking effect: for instance, the increase of upper-level clouds greatly reduces the sensitivity of the TOA fluxes to the lower-level clouds.
In summary, the non-cloud univariate feedbacks in general can be approximated well by the kernel method, although one should be mindful about the biases introduced by kernel discrepancies. One most noticeable univariate nonlinear effect in the Arctic is the surface albedo feedback. We further analyze this and other nonlinear effects in the following subsections.

Univariate Nonlinear Effects
Because radiative transfer is a complex nonlinear process (e.g., see Goody and Yung 1989), atmospheric radiation fluxes generally have a nonlinear dependency on the feedback variables and thus the univariate nonlinear effects generally exist.

Temperature
With regard to the univariate temperature feedback, a wellrecognized cause of the nonlinearity is the Planck function, although this nonlinearity is weak at the terrestrial temperatures. Based on the Stefan-Boltzmann Law, R σt 4 (9) Given the constant σ 5.67 × 10 −8 W m −2 K −4 , one may find that the nonlinear effect z 2 R zt 2 is only about 1% of the linear effect zR zt , for a 1-K perturbation around the equivalent blackbody temperature of Earth (t 255K). Another cause of the nonlinearity is the dependence of the gas absorptivity on the temperature, which also has a minor impact (Huang et al., 2007). This explains why the temperature feedback in the case of large perturbations can still be very well approximated by the linear kernels (not shown).

Water Vapor
The univariate water vapor feedback is generally well estimated when the logarithmic scaling scheme is used, although Figure 2C shows that the bias (i.e., the univariate nonlinear effect, ΔR qq ) can be non-negligible. A notable reason that causes the feedback to deviate from the logarithmic behavior is the unsaturated atmospheric absorption in the mid-infrared window around 10 μm wavelength. Here, the surface emission strongly contributes to the OLR and thus the water vapor feedback cannot be interpreted simply as the elevation of the atmospheric emission level, which gives rise to the logarithmic dependence (Bani Shahabadi and Huang., 2014). Figure 5 suggests that the water vapor feedback estimation may be improved if different scaling schemes are used for different spectral bands: logarithmic in the absorption bands (where atmospheric optical depth is large) and linear in the window bands (where optical depth is small): where the logarithmic kernel K log q accounts for the logarithmic response of the radiation flux in the absorption bands to water FIGURE 4 | Clear-sky water vapor feedback bias due to kernel bias. (A) Like Figure 2B, but using the kernels computed from different atmospheric profiles (Huang et al., 2017); (B) Bias compared to Figure 2B.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 693779 vapor perturbation and the linear kernel K lin q accounts for the linear response in the window band. Further research is required to develop and validate a global hybrid kernel set.

Surface Albedo
The univariate surface albedo feedback shows especially strong nonlinear dependence on the surface albedo a (Figure 6). This is because the multiple scattering of radiation between the surface and atmosphere renders a nonlinear dependency of planetary albedo on surface albedo. Following Stephens et al. (2015), the planetary albedo a p can be expressed as a p r + aτ 2 1 − ra (11) Here r and τ denote atmospheric reflectance and transmittance respectively; they are related as τ + r + ε 1, where ε denotes atmospheric absorptivity. This relation means that a p and thus the net shortwave radiation flux at TOA has a nonlinear dependency on a: That the radiative sensitivity to surface albedo continuously varies with the albedo value makes it difficult for any linear methods such as the kernel method to accurately measure the albedo feedback. It is interesting to notice from Eq. (12) that the radiative sensitivity decreases with a. This means that if the surface albedo kernel is computed with relatively larger albedo values under the unperturbed climate (1xCO 2 ), it will overestimate the univariate albedo feedback in a warming scenario (4xCO 2 ). This is clearly seen from Figure 6. If the kernel method is used to estimate the feedback when sea ice completely melts, the intercepts on y-axis indicate the overestimate can be serveral dozens of W m −2 . This overestimation issue was also noted in the previous studies (e.g., Block and Mauristen 2013;Zhu et al., 2019). It is also interesting to notice that although the analytical model qualitatively captures the change of radiative sensitivity to albedo, it does not accurately predict it. The neural network method proposed by Zhu et al. (2019) and Shakirova and Huang (2021) may be better suited for

Multivariate Nonlinear Effects
Besides the univariate nonlinear effects, Eq. (1) indicates that multivariate nonlinear effects, represented by such terms as zxzy ΔxΔy, may also strongly contribute to the radiative flux variations. Such terms are often referred to as the coupling effects because they result from concerted variations of the involved variables; otherwise, if their covariance were small, the average of this term over time or region would be negligible. In reality, this necessary condition is usually met because the variations of the feedback variables of concern tend to be strongly correlated. For instance, temperature warming and sea ice melt may expose more open water, which in turn leads to more evaporation, atmospheric humidity and cloudiness.
Due to large computational expenses of the brute-force RTMbased feedback calculation, we base our discussions on two representative months in the 4xCO 2 experiment: January (winter) for the longwave feedbacks and June (summer) for the shortwave feedbacks, because the two types of feedbacks are the most prominent in the two respective seasons. Figures 7,8 show the radiative feedbacks corresponding to the changes in the feedback variables illustrated in Figure 1. Table 1 summarizes their Arctic mean values. These results disclose two strongest multivariate feedback effects in the Arctic: the coupling effect between temperature and cloud, ΔR tc , in the longwave and the coupling effect between albedo and cloud, ΔR ac , in the shortwave.
In the longwave, we find that the multivariate (bivariate) feedback is dominated by the coupling effect between temperature and cloud, ΔR tc . The pattern of this coupling effect resembles, but strongly offsets, the univariate temperature feedback, ΔR (t) , which is the dominant feedback that controls the overall LW feedback in the Arctic. This coupling effect can be explained by a simple analytical model. Consider a single-layer atmosphere, with temperature t a and emissivity (absorptivity) ε, and assume the surface to be a blackbody with temperature t s : Hence, the coupling effect ΔOLR tc is found to be Because the warming in the Arctic is capped in near-surface layers, Δt a < < Δt s . This leads to reduction in OLR, offseting the increase of OLR by temperature warming. The coupling effect between temperature and water vapor can be understood in the same way. Becaue water vapor affects OLR also by affecting the atmospheric emissivity; the coupling effect ΔOLR tq is thus also affected by the factor (t 3 a Δt a − t 3 s Δt s ) in Eq. (14). Although this nonlinear effect arises from the nonuniform vertical structure of temperatuure warming and thus share the physical cause of the temperature lapse rate feedback, this nonlinear effect should be distinguished from the lapse rate feedback, which is part of the univariate temperature feedback, ΔR t . Note that the sign of ΔR tc and ΔR cq is positive in Figure 7 because the fluxes are defined to be downward positive.
In the shortwave, the dominant multivariate effect is found to be the albedo-cloud coupling, which offsets the univariate albedo feedback. From the simple model described above (Eq. 11 and Eq. 12, this can be understood as cloud-caused reduction in the atmospheric transmittance and thus reduction in the radiative senstivity to surface albedo. It is interesting to note that the patterns of some coupling feedback effects are correlated with the change patterns of the associated feedback variables. For example, the LW cloudtemperature coupling effect is correlated with surface temperature change, with a correlation coefficient of 0.86; the LW cloud-water vapor coupling effect is correlated with total water vapor (TCWV) at 0.62; the SW albedo-cloud coupling effect is correlated with the surface albedo change at 0.89. Such relation suggests that it may be possible to estimate these nonlinear effects using analytical or statistical models. Huang and Huang (2021) used such an model to explain the cloud-water vapor coupling effect and found it to be the dominant multivariate longwave feedback effect in the tropics. Although this coupling effect is not as strong Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 693779 compared to the temperature-related coupling effects in the Arctic, we find that adopting the same estimation method of Huang and Huang, 2021, Eq. 19 and using a parameter value appropriate to the Arctic (A 0.04 kg −1 m 2 ), we can very well predict the cloud-water vapor coupling effect (spatial correlation 0.99, RMSE 0.41 W m −2 ). Future works are warranted to identify methods for explaining and predicting the other coupling effects.
Lastly, for comparison, we include in Table 1 the respective feedbacks analyzed from the kernel method, i.e., the univariate linear effects for non-cloud feedbacks and the cloud feedback (ΔR p c ) obtained as a residual of total radiation change decomposition. Besides the biases in the univariate feedbacks as noted above, it is worth noting that the kernel-based estimations may also greatly bias the cloud feedback, mainly due to its negligence of the coupling effects.

CONCLUSIONS AND DISCUSSIONS
In this paper, we present an overview of the nonlinear effects in both longwave and shortwave radiative feedbacks in the CO 2 -driven Arctic warming. Based on brute-force radiation model calculations we disclose the most prominent nonlinear feedback effects and based on simple analytical models we offer explanations of their physcial causes. Although the presentation and discussion are focused on the Arctic feedbacks, the diagnostic framework (Eq. 1) and the theoretical explanations are applicable to global feedback analyses.
We identify these important nonlinear feedback effects: 1) The univariate nonlinear effect in the surface albedo feedback in the shortwave. This nonlinearity can be understood from a simple analytical model [Eq. (11)] that accounts for the coupling, due to multiple-scattering, between the surface and atmosphere (clouds). This coupling makes the radiative sensitivity to surface albedo decrease with the surface albedo value (Figure 6). Because of this effect, it generally leads to an overestimate of the surface albedo feedback when albedo kernels computed from the current climate are used to quantify the albedo feedback in a warming climate ( Figure 3; Table 1).
2) The bivariate surface albedo-cloud coupling effect in the shortwave. This effect is attributable to the masking effect of cloud increase that damps the radiative sensitivity to surface albedo. This effect is the most prominent in the summer when solar insolation is strong, as illustrated by Figure 8 for the month of June.
3) The multivariate temperature-cloud feedback in the longwave. This effect is attributable to the fact that the Arctic warming is much stronger at and near the surface than in the upper air, which leads to a damping effect on the temperature feedback. This nonlinear effect should be distinguished from the temperature lapse rate feedback and is found to be the strongest in the winter as illustrated by Figure 7 for the month of January. 4) Although the univariate water vapor feedback largely scales logarithmically with water vapor changes, it is found that the relation deviates from the logarithmic scaling, especially in the window band. This is due to the unsaturated atmospheric absorption in this band and suggests that a hybrid scaling method as proposed by Eq. (10) may improve the accuracy of the kernel-diagnosed water vapor feedback (compare Figure 2C and Figure 5C).
It should be noted that the large nonlinear effects discovered here is not limited to the 4xCO 2 experiment. As shown by Huang and Huang (2021) for the longwave feedbacks and Shakirova and Huang (2021) for the shortwave feedbacks, similar, strong nonlinear effects exist even also interannual climate variations. It is noted that the nonlinear effects may quantitatively differ in different forcing experiments, thus requiring them to be assessed more comprehensively in future work.
The strong nonlinear effects as disclosed here call into question the accuracy of linear methods currently used in the feedback analysis. Nonlinear methods are needed to improve the accuracy of feedback quantificaiton when RTM-based PRP experiments are not feasible due to its forbidding computational demands. Especially in need are replacement of the linear kernels for the surface albedo feedback and cloud feedback quantification. Although a handful of TABLE 1 | Arctic mean all-sky feedbacks in the 4xCO2 experiment for the two selected months. Units: W m-2. Area-weighted averages are taken for the region 70-90°N. Two sets of radiative kernels have been used to measure the univariate linear feedbacks: Ker1 is computed from the GCM instantaneous profiles in this work and thus is of no kernel bias; Ker2 is the kernel computed by Huang et al. (2017) from the ERA-interim reanalysis profiles, which leads to biases in diagnosed univariate feedback as explained by Eq. 8. LW (jan ΔR (t,q,c) ΔR ( studies have touched this topic, for instance, using quadratic fitting , histogram (Zelinka et al., 2012) and neural network (Zhu et al., 2019) methods, this challenging problem demands devoted research programs to further develop, test and mature the candidate methods.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: The CESM codes can be downloaded from National Center for Atmospheric Research (NCAR) website (http://www.cesm.ucar.edu/models/cesm1.2/). The RRTM code can be downloaded at http://rtweb.aer.com/rrtm_frame.html.

AUTHOR CONTRIBUTIONS
YH designed the research and wrote the paper. HH conducted the longwave feedback analysis and AS conducted the shortwave feedback analysis.