Effects of Anthropogenic and Natural Forcings on the Summer Temperature Variations in East Asia during the 20th Century

The effects of the emissions of anthropogenic greenhouse gases (GHGs), aerosols, and natural forcing on the summer-mean surface air temperature (TAS) in the East Asia (EA) land surface in the 20th century are analyzed using six-member coupled model inter-comparison project 5 (CMIP5) general circulation model (GCM) ensembles from five single-forcing simulations. The simulation with the observed GHG concentrations and aerosol emissions reproduces well the land-mean EA TAS trend characterized by warming periods in the early (1911–1940; P1) and late (1971–2000; P3) 20th century separated by a cooling period (1941–1970; P2). The warming in P1 is mainly due to the natural variability related to GHG increases and the long-term recovery from volcanic activities in late-19th/early-20th century. The cooling in P2 occurs as the combined cooling by anthropogenic aerosols and increased volcanic eruptions in the 1960s exceeds the warming by the GHG increases and the nonlinear interaction term. In P3, the combined warming by GHGs and the interaction term exceeds the cooling by anthropogenic aerosols to result in the warming. The SW forcing is not driving the TAS increase in P1/P3 as the shortwave (SW) forcing is heavily affected by the increased cloudiness and the longwave (LW) forcing dominates the SW forcing. The LW forcing to TAS cannot be separated from the LW response to TAS, preventing further analyses. The interaction among these forcing affects TAS via largely modifying the atmospheric water cycle, especially in P2 and P3. Key forcing terms on TAS such as the temperature advection related to large-scale circulation changes cannot be analyzed due to the lack of model data.


Introduction
Intense research in the recent decades confirmed that the increase in the global-mean surface air temperature (TAS), most notably since early 20th century, has been caused by the emissions of anthropogenic greenhouse gases (GHGs) and aerosols/aerosol precursors from industrial sources with minimal contributions from the natural variability [1][2][3][4][5][6]. The human-induced climate change affects a wide range of sectors, both human and natural, over the globe. Because of wide-spread impacts on the human society and the natural ecosystems, understanding the ongoing climate variations and projecting future climate have become an important practical concern worldwide [7], especially for to more realistic variations in the anthropogenic GHG, aerosols, natural forcings, and the internal variability of the climate system, using transient simulations. Song et al. [35] attempted to study the effects of the observed variations in the GHGs and aerosol forcing on the EA climate using CMIP5 simulations for three different simulations with historical, historical GHG, and historical natural forcing; however, their study based on the linearity assumption cannot separate the nonlinear interactions from the aerosol forcing. This study is designed to explore the effects of GHGs, anthropogenic aerosols, and natural forcing explicitly with separate aerosols forcing runs, so that the effects of interactions among these forcing in shaping the observed TAS variations over EA in the 20th century can be account for. We have analyzed a total of 30 simulations from 6 CMIP5 models for five different simulations with historical, historical GHGs, historical anthropogenic aerosols, and historical natural forcings, and a pre-industrial (PI) conditions in order to assess separately the effects of anthropogenic GHGs and aerosols, as well as natural forcings and internal variability on TAS over EA. The CMIP5 models and experiments setup are described in Section 2. It is followed by the analyses of the model data in Section 3. Summaries and Discussions are presented in Section 4.

Model Data
In this study, we use 30 simulations from 6 GCMs contributing to phase 5 of the Coupled Model Intercomparison Project (CMIP5) through the Earth System Grid Federation following Taylor et al. [36]. The model names, institutions, and resolutions are listed in Table 1. All the model outputs are bilinearly interpolated into horizontal resolution of 2.81 • × 2.81 • grid, which is the same as grid of the Canadian Earth System Model 2 (CanESM2) model. The bilinear interpolation tends to underestimate local maxima but does not create fictitious local peaks, either. The multi-model ensemble (MME) is examined by the arithmetic mean of model output, with the same weight for each model. Five experiments of different forcing are selected for each model: the historical simulation (HST) is forced by observed history of the anthropogenic GHGs and aerosols as well as the natural forcing by the solar and volcanic activities. The single-forcing runs using only the historical GHGs (GHG), anthropogenic aerosol emissions (AER), and natural forcing (NAT) are used to estimate the effects of these individual forcings on TAS. The effects of the internal variability of the climate system (INT) is estimated from the Pre-industrial (PI) simulation in which the solar constant is fixed as a constant, no volcanic activities are considered, and the GHG and aerosol concentration are fixed at the PI level. Each GCM provides a single member to each ensemble. More details are described below Section 2.3. GFDL-ESM2M [40] Same as GFDL-CM3 (NOAA/GFDL) 2.5 • × 2 • IPSL-CM5A-LR [41] Institute Pierre-Simon Laplace (IPSL) 3.75 • × 1.87 • NorESM1-M [42] Norwegian Climate Centre (NCC) 2.5 • × 1.875 •

Observational Data
This study employs temperature indices derived from high-resolution (0.5 • × 0.5 • ) monthly-mean temperatures from the Climate Research Unit (CRU) TS version 3.2.2. [43] for the reference for evaluating the simulated land-mean temperatures. The CRU dataset is based on analyses of temperature data at surface stations and covers the global land surface except Antarctica. The CMIP5 simulations used in this study account the temporal variations in the solar and volcanic activities towards the natural forcing. The solar activities are prescribed from the annual mean time series of sunspot number from the Solar Influences Data Analysis (SIDC, Royal Observatory of Belgium) and the annual mean time series of total solar irradiance data from NOAA [44][45][46]. The effects of volcanic eruptions on TAS are accounted in terms of the radiative effects of volcanic ashes by prescribing the monthly mean time series of the stratosphere aerosol optical depth (SAOD) dataset from the National Aeronautics and Space Administration/Goddard Institute for Space Studies (NASA/GISS) [33,47,48].

Methodology
This study explores the effects of anthropogenic (by the emissions of aerosols and GHGs) and natural (by the solar and volcanic activity) forcings on the historical variations of the summer TAS in EA in the 20th century. For a given climate variable Ψ, the amount of variations over a time period (∆Ψ) that are associated with the changes in anthropogenic aerosol concentration (∆α), the GHGs concentrations (∆g), and the natural forcings (∆n) may be decomposed as in Equation (1) below: ∆Ψ(∆α, ∆g, ∆n) = (∂Ψ/∂α) ∆a + (∂Ψ/∂g) ∆g + (∂Ψ/∂n) ∆n + R(∆α, ∆g, ∆n) + I where R(∆α, ∆g, ∆n) is the change due to the nonlinear interactions between ∆α, ∆g and ∆n, and I is the changes due to internal variability. The first three terms on the right hand side of Equation (1) represent the changes in Ψ solely by the forcing from anthropogenic aerosols, GHGs, and natural events (the solar and volcanic activities) in order. The five term in Equation (1) are estimated using the six-model ensembles from five sensitivity simulations for the period 1860-2005 ( Table 2). The first simulation from which ∆Ψ(∆α, ∆g, ∆n) is calculated accounts for the historical variations in the emissions and concentrations of aerosols, GHGs, and natural forcing; i.e., the most realistic representations of all forcings (HST) on TAS. Details of the forcing agents (GHGs, Aerosols, and Natural, etc.) in CMIP5 models are described in Table S1 of the supporting information. Three single-forcing simulations ( Table 2) are used to calculate the first three terms on the right hand side of Equation (1). The anthropogenic aerosol-only simulation (AER) accounts for the historical anthropogenic aerosol emissions while the GHGs concentrations are fixed at the 1860 level. The GHG-only simulation (GHG) implements the historical GHGs concentration while anthropogenic aerosol emissions are fixed at the 1860 level. The natural forcing run (NAT) implements the historical natural forcing due to the solar activity and volcanic eruptions while the GHGs and aerosol concentrations are fixed at the 1860 level. In addition, the PI simulation (INT), where both the concentration of the anthropogenic aerosols and GHGs as well as the natural forcings are fixed at the pre-industrial (1860) level, is employed to assess the effects of the internal variability of the climate system on the TAS variations.
Because internal variabilities of various time scales can play an important role in modulating regional temperature and characterizing the individual model uncertainties, we first examined the significance of internal variabilities of 30-year and 150-year timescales using the Pi-control runs. The trend of the 6-model ensemble temperature over the 150-year period is much smaller than those from the simulations listed in Table 2 and the corresponding residual term (not shown). For estimating the trends at 30-year time scales, we randomly selected 100 samples of consecutive 30-year periods and calculated 30-year trends for individual samples. The 30-year trends calculated for the 100 samples are used to calculate the mean and the standard deviation of the 30-year trend. This step was necessary because the Pi-control runs are not tied to real calendar, i.e., the Pi-control runs cannot be related to specific periods. The mean 30-year trends over the 100 samples are nearly zero; i.e., negligible compared to the effects obtained in the simulations and the corresponding residuals. The uncertainty of the 30-year trends due to the internal variability, estimated by the standard deviation of the 30-year trends of the 100 samples is also smaller than most terms except the effects of aerosols in P1 and the effects of the natural forcing term in P3. The uncertainty range is presented by the two dashed black lines in Figure 1c. Note that the effects of aerosols and the natural forcing in P1 and P3, respectively, are negligible compared to the effects of other forcings in the same period. Based on the analysis above, we concluded that the effects of internal variability on the temperature trend in each period can be neglected. Table 2. Five sets of climate simulations are used in this study.

HST
The twentieth century historical climate simulations (historical; all forcing runs). The forcing data are from the CMIP5 for IPCC AR5. The scenario data include natural agent change (mainly solar and volcanic) and anthropogenic agent changes (well-mixed GHGs, aerosols, ozone, land-use changes)

AER
Aerosols-only forcing simulations (historicalMisc runs). Note that the level of complexity and the range of aerosol species differ among different CMIP5 models (see Table S1 of the supporting information)

GHG
Greenhouse gases-only forcing simulation (historicalGHG runs). The simulations are forced under well-mixed GHGs as in the all forcing run, but other agents are fixed at the pre-industrial level NAT Natural-only forcing simulation (historicalNat runs). The natural agents are the same as those used in the all forcing run, but anthropogenic agents are fixed at the pre-industrial level

INT
Pre-industrial simulation (piControl runs). Both anthropogenic and natural forcings are fixed at the pre-industrial level A detailed description of the CMIP5 simulation design is provided by Taylor et al. [36]. These five simulations allow us to estimate all five terms in Equation (1). The residual term R (the fourth term in the right-hand-side of Equation (1)) that represents the nonlinear interaction of anthropogenic (GHG, aerosol) and natural forcings is calculated by subtracting the sum of the single forcing experiments AER, GHG, NAT and INT from HST. The CMIP5 model data are analyzed for the East Asia region covering 20 • -50 • N and 110 • -140 • E which includes eastern China, Korea, Japan, and Manchuria. The summer climate is represented by the average over the three summer months, June, July, and August. Statistical significance level is determined based on the Mann-Kendall trend test. Based on the observed TAS variations in Figure 1b, we select three 30-year periods based on the TAS trend variations: 1911-1940 (P1: warming), 1941-1970 (P2: cooling) and 1971-2000 (P3: warming). We explore the effects of the anthropogenic GHGs and aerosols, natural forcings, as well as their interactions, on the temperature records for each period from the simulated data. The variations of the selected variables in the three periods are represented by the trends over the corresponding periods.

Evaluation of the 20th Century Overland Temperature in HST
The interannual variations of the land-mean TAS in HST are compared against the CRU data for the entire globe ( Figure 1a) and for EA (Figure 1b). HST is run with the observed emissions of well-mixed GHGs, anthropogenic aerosols, and natural agent change (solar and volcanic) during the simulation period; thus is expected to generate the most realistic TAS among the five simulations. The observed TAS over the EA land region (Figure 1b) is in-phase with the global TAS (Figure 1a), but with more pronounced variability during the middle  of the 20th century period. For both the global and EA land surfaces, HST well depicts the long-term TAS variations characterized by a positive TAS anomaly peak in the 1940s and a negative peak in early 1970s followed by continuous increases in the late 20th century period. The model performs better for the global mean than for the EA-mean; the correlation coefficient between the observed and the TAS time series in HST is 0.95 for the global land surface and 0.68 for the EA land surface. The five sensitivity simulations show that the simulated TAS variations during the 20th century primarily result from the warming effects of GHG increases and the cooling effects of anthropogenic aerosols (red and blue lines, respectively, in Figure 1a,b). The natural forcing by the solar activity and volcanic eruptions also contributed to the warming trend in early 20th century and the cooling trend in mid-20th century (green lines in Figure 1a,b). Details of the TAS trends in individual sensitivity experiments will be discussed in the following sub-sections.
Atmosphere 2019, 10, x FOR PEER REVIEW 6 of 23 (green lines in Figure 1a,b). Details of the TAS trends in individual sensitivity experiments will be discussed in the following sub-sections.  Table 2 and the corresponding residuals for the three periods. In (a) and (b), the 10 year running-mean TAS anomaly  over the entire globe and EA are shown. CRU indicates the CRU data and RES is the residual term calculated from the simulations as Equation (2); abbreviations for other runs are the same as in Table 2 Table 2 and the corresponding residuals for the three periods. In (a) and (b), the 10 year running-mean TAS anomaly  over the entire globe and EA are shown. CRU indicates the CRU data and RES is the residual term calculated from the simulations as Equation (2); abbreviations for other runs are the same as in Table 2. Error bars indicate the standard deviation of CMIP5-model multi-model ensemble (MMEs) from each simulation. The black dashed horizontal lines indicate the range of the internal variability on the 30 year timescale estimated as one standard deviation of the 30-year trends of the 100 samples from the Pi-Control MME as described in Section 2.

TAS Trends Related to Anthropogenic Emissions
The TAS trend in the four sensitivity simulations (HST, AER, GHG, and NAT) are compared against the observed trend for the three periods defined in Section 3.1 (Figure 1c). The mean 30-year internal variability estimated from the Pi-control runs is negligible; the range of the TAS trend due to internal variability is presented by the black dashed lines in Figure 1c. The TAS trend in HST is different from the sum of the trends in the single-forcing simulations. These differences arise due to the interactions among the individual forcing elements that cannot be simulated in the single-forcing runs. The residual term R in Equation (1) that arises from interactions among the individual forcing elements, is estimated from the five sensitivity simulations as in Equation (2) below: Note that the internal variability term (INT) is neglected in calculating the RES term in (2) as its contribution to the TAS trend is negligible compared to the contribution from other forcings. The sign of the linear trend of TAS in HST (purple bar in Figure 1c) agrees with that derived from the CRU data (black hatched bar in Figure 1c) for all of the three 30-year periods despite some differences in the magnitude. Results from the historical (HST) and the three single-forcing simulations (GHG, AER, and NAT) and the residual term (RES) obtained from Equation (2) are analyzed to assess the effects of GHGs, anthropogenic aerosols, natural forcings and their nonlinear aspects on the TAS variations during the 20th century in the following sub-sections.

TAS Variations in P1
The first 30-year period 1911-1940 (P1) is characterized by a rapid increase in TAS (Figure 1b), mainly due to the natural forcing (increasing solar irradiance and decreasing volcanic aerosols) and the increase in GHGs as the main contributors to the positive TAS trend (Figure 1c). Anthropogenic aerosols exert almost negligible effects on TAS in P1. The interaction term (RES) is negative in P1 but is smaller in the magnitude and is of much larger statistical uncertainty compared to GHG and NAT. Thus, the positive TAS trend during P1 mainly results from the combined warming effects of the natural forcing and GHG increases exceeding the cooling effects of the interaction term. The global CO 2 concentration (red dashed line) and EA SO 2 emission (blue solid line) profiles implemented in the HST, AER and GHG runs show that both the CO 2 concentration and the increase in the East Asian SO 2 emission is in P1, a pre-industrial period for most of EA, is much lower than in P2 and P3 in which heavy industrialization has occurred in EA ( Figure 2). This results in smaller effects of GHGs and aerosols on the TAS variations in P1 (early 20th century) than in P2/P3 (middle/late 20th century). The positive contribution of the natural variability to the warming trends in P1 is a part of the long-term recovery of TAS from the large cooling by the massive volcanic activities in late 19th/early 20th century (the Krakatoa 1883, the Santa Maria 1902, the Mt. Katmai 1912 in Figure 3b) [45,47,49]. The solar activity has also increased in P1, but the magnitude is small. Details of the solar activity and volcanic eruptions implemented in the CMIP5 natural forcing are referred to Moss et al. [50].     Warming trends in P1 are observed in most of the EA land region, particularly in the Manchuria, the Korean Peninsula, far-eastern Russia, and Japan, except the Mongolia region (Figure 4a). The simulated TAS trends in P1 show large regional variations over EA. The simulated TAS in HST also shows positive trends over the land in P1 but the trends are smaller than from the CRU data ( Figure 4b). These TAS trend closely coincides with the TAS response to the natural forcing (NAT), especially over northeastern and southern China (Figure 4e). For the East Sea of Korea, much of the warming by GHG increases (Figure 4d) is offset by the cooling by anthropogenic aerosols (Figure 4c). Note that statistically significant TAS trend in both AER and GHG occurred mainly over the seas around Japan (Figure 4c,d). The interactions between the anthropogenic and natural forcings (Figure 4f) result in very weak negative (positive) TAS trends over land (ocean).
Atmosphere 2019, 10, x FOR PEER REVIEW 9 of 23 Warming trends in P1 are observed in most of the EA land region, particularly in the Manchuria, the Korean Peninsula, far-eastern Russia, and Japan, except the Mongolia region (Figure 4a). The simulated TAS trends in P1 show large regional variations over EA. The simulated TAS in HST also shows positive trends over the land in P1 but the trends are smaller than from the CRU data ( Figure  4b). These TAS trend closely coincides with the TAS response to the natural forcing (NAT), especially over northeastern and southern China (Figure 4e). For the East Sea of Korea, much of the warming by GHG increases (Figure 4d) is offset by the cooling by anthropogenic aerosols (Figure 4c). Note that statistically significant TAS trend in both AER and GHG occurred mainly over the seas around Japan (Figure 4c-d). The interactions between the anthropogenic and natural forcings (Figure 4f) result in very weak negative (positive) TAS trends over land (ocean). The surface downwelling radiation budget is analyzed to elucidate the role of radiative forcing on TAS. It must be noted that TAS is shaped not only by radiative forcing but also feedback within the climate system that involves circulation changes that accompany temperature advection. Despite their importance, the contribution from the climate feedback cannot be analyzed in this study due to the lack of suitable data. In the early 20th century period (P1), changes in the total aerosol optical thickness (contour line in Figure 5b) is small, resulting in weak SW trends in the clear-sky condition (Figure 5b, shading). This is consistent with the result that the TAS trend by anthropogenic aerosols The surface downwelling radiation budget is analyzed to elucidate the role of radiative forcing on TAS. It must be noted that TAS is shaped not only by radiative forcing but also feedback within the climate system that involves circulation changes that accompany temperature advection. Despite their importance, the contribution from the climate feedback cannot be analyzed in this study due to the lack of suitable data. In the early 20th century period (P1), changes in the total aerosol optical thickness (contour line in Figure 5b) is small, resulting in weak SW trends in the clear-sky condition ( Figure 5b, shading). This is consistent with the result that the TAS trend by anthropogenic aerosols is almost negligible in P1. Much of the SW radiation trends (shading in Figure 5c) occur under cloudy sky conditions in close negative correlation with the changes in cloudiness (contour line in Figure 5c). The land-mean SW trend in the four runs show that the positive SW trend by the natural variability (NAT) is over ridden by the negative SW trend by anthropogenic aerosols (AER) and greenhouse gases (GHG) in Table 3. In conjunction with the increase in cloudiness (Table 4), the SW trend in HST becomes negative, i.e., the SW budget works to cool TAS. The cooling effects of the SW budget is compensated by large positive effects of the LW budget to contribute to the positive TAS trend in HST during P1. It is difficult to disseminate the LW trend in Table 3 into the one due to the forcing by each element and the one due to the feedback through the climate system. However, the large LW trend in GHG in the clear-sky condition suggests that much of the positive LW trend in HST can be assumed as an indicator of the atmospheric greenhouse effect as experienced at the surface [52][53][54][55][56], and thus a forcing towards the positive TAS trend in P1. Table 3 also shows that the SW trend in the cloudy-sky condition exceeds that in the clear-sky condition. In conjunction with the large positive cloudiness trend in P1 (Table 4), the negative SW trend in HST is strongly affected by the increase in cloudiness.
Atmosphere 2019, 10, x FOR PEER REVIEW 10 of 23 is almost negligible in P1. Much of the SW radiation trends (shading in Figure 5c) occur under cloudy sky conditions in close negative correlation with the changes in cloudiness (contour line in Figure 5c). The land-mean SW trend in the four runs show that the positive SW trend by the natural variability (NAT) is over ridden by the negative SW trend by anthropogenic aerosols (AER) and greenhouse gases (GHG) in Table 3. In conjunction with the increase in cloudiness (Table 4), the SW trend in HST becomes negative, i.e., the SW budget works to cool TAS. The cooling effects of the SW budget is compensated by large positive effects of the LW budget to contribute to the positive TAS trend in HST during P1. It is difficult to disseminate the LW trend in Table 3 into the one due to the forcing by each element and the one due to the feedback through the climate system. However, the large LW trend in GHG in the clear-sky condition suggests that much of the positive LW trend in HST can be assumed as an indicator of the atmospheric greenhouse effect as experienced at the surface [52][53][54][55][56], and thus a forcing towards the positive TAS trend in P1. Table 3 also shows that the SW trend in the cloudy-sky condition exceeds that in the clear-sky condition. In conjunction with the large positive cloudiness trend in P1 (Table 4), the negative SW trend in HST is strongly affected by the increase in cloudiness.

TAS Variations in P2
The temperature variations over EA in the second period (P2) are characterized by a notable cooling trend (Figure 1c). The linear trends of TAS in each simulation show that the observed cooling in EA during P2 is a consequence of the combined cooling effects of anthropogenic aerosols and the natural forcing that exceeding the combined warming effects of the GHGs increases and the interactions (Figure 1c). The strong negative TAS trend in the CRU data occurs in southern and northeastern China, and Mongolia, with positive trends in eastern China, Taiwan and Japan (Figure 6a). The TAS trend in HST is weaker than the TAS trend from the CRU data (Figure 1c, Figure 6b). These negative TAS trend closely coincides with the TAS response to the anthropogenic aerosols and the natural forcing (Figure 6c (Figure 6d). The positive TAS trends over the EA land surface, mostly due to the GHG increase and the interactions, offset the cooling trends in south China, Mongolia, the coastal areas around China and Korea, and the ocean between Japan and far-eastern Russia.
The regional variations of the TAS trend in HST (Figure 6b) are closely related to SW variations ( Figure 7a) with a spatial pattern correlation of 0.79 between the TAS trend and the SW trend. The all-sky SW trend is determined mainly in the clear-sky condition (shading in Figure 7b) and is closely correlated with the total aerosol optical depth trend (contour line in Figure 7b); however, unlike in P1, SW reduction occurs more clearly in the clear-sky condition than in the cloudy-sky condition; i.e., the direct radiative forcing by the rapidly increasing SO 2 emissions and the associated increase in sulfate aerosols has played an important role in the overall cooling trend during P2. The positive SW trend in the ocean area (shading in Figure 7a,c) is due to the decrease in total cloudiness (contour line in Figure 7c). The LW trend in HST is negative over the EA region (Figure 7d), mostly under the clear sky condition (Figure 7e). This is somewhat peculiar because the GHGs concentration increases rapidly ( Figure 2) and works to increase TAS (Figure 1c; Figure 7c) during P2. This will be further examined below using the land-mean LW trend in individual single-forcing simulations.  Unlike for P1 (and P3 to be discussed in the following section), the negative land-mean SW trend in HST is consistent with the TAS trend in this period. The land-mean SW and LW trends from each single forcing runs for this period (Table 5) show that AER and the NAT forcing are related to large negative SW trend while the GHG forcing provides large positive LW trend while the LW trend in AER and NAT is negative. Thus, the negative LW trend in HST largely results from the combined negative LW effects in AER and NAT dominating the positive effects of in GHG. Similar to the SW trend, the LW trend in P2 are dominated in the clear-sky conditions in all runs except the SW trend by NAT. Table 4 shows that the magnitude of the land-mean cloudiness trend in HST is smallest in P2 among the three periods and that the natural forcing yields large positive cloudiness trend. Unlike for P1 (and P3 to be discussed in the following section), the negative land-mean SW trend in HST is consistent with the TAS trend in this period. The land-mean SW and LW trends from each single forcing runs for this period (Table 5) show that AER and the NAT forcing are related to large negative SW trend while the GHG forcing provides large positive LW trend while the LW trend in AER and NAT is negative. Thus, the negative LW trend in HST largely results from the combined negative LW effects in AER and NAT dominating the positive effects of in GHG. Similar to the SW trend, the LW trend in P2 are dominated in the clear-sky conditions in all runs except the SW trend by NAT. Table 4 shows that the magnitude of the land-mean cloudiness trend in HST is smallest in P2 among the three periods and that the natural forcing yields large positive cloudiness trend.

TAS Variations in P3
For P3, HST shows positive TAS trend of a similarly in magnitude as in the CRU data over the EA land surface (Figure 1c). The TAS trend in GHG continues to increase (red line in Figure 1b) in response to the the rapid increase of carbon dioxide concentration (Figure 2) in the second half of the 20th century. Unlike in GHG, the negative TAS trend in AER has slowed since the mid-1980s (blue line in Figure 1b,c). Nevertheless, the negative TAS trends by anthropogenic aerosols is the only measurable offset for the positive TAS trend due to the GHGs effect (Figure 1c). Natural forcing also exerts negative TAS trend; however, its effects on the TAS trend in P3 are within the range of the uncertainty associated with the internal variability, and are insignificant compared to the effects of GHGs and anthropogenic aerosols. The weak effects of natural forcing are due to the absence of strong volcanic eruption (except Pinatubo 1991) during late 20th century (Figure 3b). It is found that much of the cooling trend due to natural forcing has been recovered since 1991 in the TAS time series for EA as well as for the entire globe (green line in Figure 1a,b). This result is consistent with the satellite-based estimates of solar irradiance changes (−0.1 to +0.1 Wm −2 ) associated with volcanic eruption from 1980 to 2010 [5].
The warming trend in the CRU data occurs in most of the EA region, particularly in northern EA including Mongolia, the Manchuria, the Korean Peninsula, and Japan (Figure 8a), with cooling trends in central China. The warming in HST also shows strong positive TAS trend of a similarly in magnitude and spatial pattern as in the CRU data (Figure 8b). The positive TAS trend is mostly due to GHGs increases (Figure 8d) with partial effects of the interaction over northeastern China, far-eastern Russia, and Western North Pacific (WNP) (Figure 8f). The negative TAS trend in central China closely coincides with the TAS response to the anthropogenic aerosols (Figure 8c). The anthropogenic aerosols exert significant decreasing TAS trends in central China, northern Korea, the South China Sea, and the WNP. The natural forcing (Figure 8e) yields weak decreasing TAS trend within the domain. The negative TAS trends from anthropogenic aerosols and natural forcing works to offset the warming trend due to the GHGs and interaction term in the late 20th century as noticed in the land-mean temperature trend (Figure 1c). Like in P2, the negative SW trend (Figure 9a) occurs largely in the clear-sky condition (shading in Figure 9b), and is closely correlated with the total AOD trend (contour line in Figure 9b). The AOD trend is positive in the domain except Japan (contour in Figure 9b), contributing to the negative clearsky SW trend along the heavily industrialized east coast region of China (shading in Figure 9b). Like in P2, the negative SW trend (Figure 9a) occurs largely in the clear-sky condition (shading in Figure 9b), and is closely correlated with the total AOD trend (contour line in Figure 9b). The AOD trend is positive in the domain except Japan (contour in Figure 9b), contributing to the negative clear-sky SW trend along the heavily industrialized east coast region of China (shading in Figure 9b). Reduced cloudiness in Mongolia, South Korea, and the Far-east Russia (contour red dotted line in Figure 9c) locally contributes to positive SW trends (Figure 9c, shading). Thus, the rapid increase in anthropogenic aerosols affect TAS mostly through aerosol direct effects in this period. The LW trend is positive (Figure 9d) and is mostly in the clear-sky condition (Figure 9e). Like in P1 and P2, the geographical variations of the TAS trend in HST closely resemble the LW trend in the clear-sky condition (Figure 9d).  Table 4 shows that the feedback through clouds in P3 may be smaller than that in P1 although it is larger than in P2.   The land-mean SW and LW trends during P3 are similar to those in P1 as shown in the trends in individual simulations (Table 6). For example, despite the positive TAS trend in HST, the SW trend is negative. The leading terms for the SW and LW trend are the negative trend in AER and the positive trend in GHG, respectively. Thus, the effects of negative SW of anthropogenic aerosols is overcome by the positive effects of GHG increases to contribute to the positive TAS trend in P3. Table 4 shows that the feedback through clouds in P3 may be smaller than that in P1 although it is larger than in P2. Table 6. The land-mean linear trend (Wm −2 decade −1 ) of the JJA-mean downwelling SW, LW, and the net radiation at the surface from the HST, AER, GHG, and NAT simulations for P3. The numbers in the parenthesis are the corresponding values for the clear-sky and cloudy-sky condition.

Effects of Interactions among Individual Forcing Elements
The residual term RES in (2) arises from the feedback among the forcings by greenhouse gases, aerosols, and natural variability that cannot be estimated in single-forcing experiments. As shown in Figure 1c, the feedback among these forcings to the climate system is also important in determining the TAS trend. According to Feichter et al. [23], the feedback through the atmospheric hydrologic cycle plays an important role in shaping the nonlinear aspect of the climate response to GHGs and aerosol forcings. Because it is practically impossible to isolate individual feedback processes in climate model simulations, we attempt to estimate the effects of the interaction term relative to the sum of the effects from the single forcing simulations. The interaction term in P1, the early 20th century period, is of smaller/larger magnitude/uncertainty compared to the sum of the effects from all single-forcing runs. Since the mid-20th century, the effects of the interaction term in the TAS trend become more obvious with much smaller uncertainty compared to that in P1 (Figure 1c). To identify plausible processes involved in the interaction term, we compare the linear trends of hydrological variables (such as water vapor, cloudiness, and precipitation rate etc.) in the all-forcing runs (HST) against those in the sum of single-forcing simulations ( Figure 10): Examinations of the water vapor path, cloudiness, cloud-water path, and precipitation show that the interaction term works to decrease the atmospheric water vapor and cloud water contents and cloudiness as well as precipitation in P2 ( Figure 10). The decreasing trend in the condensed phase of hydrometeors, cloud water (−346% of HST) and precipitation (−22% of HST), is much larger than in water vapor (−13% of HST). The reduced atmospheric water vapor, in turn, decreases condensation to enhance the negative trends in the cloud water path and precipitation (Figure 10c,d). The negative trend in the total precipitation in HST is mainly by decreased convective cloud (Figure 10e), whereas the stratiform precipitation trend is almost negligible in P2 (purple in Figure 10f). The large decrease in cloudiness and cloud water path (Figure 10b,c) by the interaction term suggests that the decreased water vapor path results in suppressing both convective precipitation (−9% of HST) and stratiform precipitation (−280% of HST) (Figure 10e,f); however, the amount of stratiform precipitation is much smaller than the amount of convective precipitation.
In P3, the interaction term contributes to the increase in atmospheric water vapor, cloud water contents, and precipitation, but exerts minimal effects on cloudiness (Figure 11a-d). The increase in the condensed phase of hydrometeors, cloud water (46% of HST) and precipitation (88% of HST), is much larger than in water vapor (24% of HST). The enhanced atmospheric water vapor, in turn, increases condensation to enhance the positive trends in cloud water path ( Figure 11c) and precipitation (Figure 11d). Like in P2, the positive trend in the total precipitation in HST is determined mainly by the convective cloud precipitation (Figure 11e); unlike in P2, the interaction term contribute to the large increase in convective precipitation (133% of HST). The slight decrease (almost negligible, −3% of HST) in cloudiness by the interaction term (Figure 11b) suggests that the increased water vapor and cloud water contents by the interaction term largely result in enhancing deep clouds than stratiform clouds (−133% of HST). (Figure 11e,f). Note that the descriptions on the effects of the interaction term here are largely qualitative. Legitimate quantitative analyses are not feasible because of the complexity in the feedbacks and also because it requires much more data that are available from the CMIP5 archive.
cloudiness as well as precipitation in P2 ( Figure 10). The decreasing trend in the condensed phase of hydrometeors, cloud water (−346% of HST) and precipitation (−22% of HST), is much larger than in water vapor (−13% of HST). The reduced atmospheric water vapor, in turn, decreases condensation to enhance the negative trends in the cloud water path and precipitation (Figure 10c,d). The negative trend in the total precipitation in HST is mainly by decreased convective cloud (Figure 10e), whereas the stratiform precipitation trend is almost negligible in P2 (purple in Figure 10f). The large decrease in cloudiness and cloud water path (Figure 10b,c) by the interaction term suggests that the decreased water vapor path results in suppressing both convective precipitation (−9% of HST) and stratiform precipitation (−280% of HST) (Figure 10e,f); however, the amount of stratiform precipitation is much smaller than the amount of convective precipitation. In P3, the interaction term contributes to the increase in atmospheric water vapor, cloud water contents, and precipitation, but exerts minimal effects on cloudiness (Figure 11a-d). The increase in the condensed phase of hydrometeors, cloud water (46% of HST) and precipitation (88% of HST), is much larger than in water vapor (24% of HST). The enhanced atmospheric water vapor, in turn, increases condensation to enhance the positive trends in cloud water path ( Figure 11c) and precipitation (Figure 11d). Like in P2, the positive trend in the total precipitation in HST is determined mainly by the convective cloud precipitation ( Figure 11e); unlike in P2, the interaction term contribute to the large increase in convective precipitation (133% of HST). The slight decrease (almost negligible, −3% of HST) in cloudiness by the interaction term ( Figure 11b) suggests that the increased water vapor and cloud water contents by the interaction term largely result in enhancing deep clouds than stratiform clouds (−133% of HST). (Figure 11e,f). Note that the descriptions on the effects of the interaction term here are largely qualitative. Legitimate quantitative analyses are not feasible because of the complexity in the feedbacks and also because it requires much more data that are available from the CMIP5 archive.

Summary and Discussion
The relationships between the multi-decadal variations of the summer-mean TAS and the emissions of anthropogenic GHGs, aerosols, and natural agents over the EA land surface in the 20th century are examined using six-member CMIP5 GCM ensembles from five sensitivity simulations in which various combinations of the historical variations of GHGs, anthropogenic aerosol emissions, and the natural variability are implemented. Unlike most previous studies that explore these effects for an equilibrium state [23,34], we attempt to understand the relationship between the historical TAS variation and these external forcing more realistically from transient simulations. We also examined the effects of the internal variability on the summer TAS over EA using the simulations in which GHG concentrations and anthropogenic aerosol emissions are fixed at the pre-20th century level (INT). For the analysis period, both observed and simulated summer TAS trends are mostly stronger than the magnitude of internal variability, suggesting that the observed summer TAS variations are primarily induced by anthropogenic and natural forcing with minimal effects from the internal variability.
The CRU data shows that the 20th century TAS variations are characterized by warming trends for the early  and late  20th century periods with a period of cooling trend between them . The six-member HST-run ensemble that includes the historical variations of GHG concentrations, anthropogenic aerosol emissions, and natural forcings by solar and volcanic activities, simulates reasonably well the observed variations of the land-mean summer TAS during the 20th century for EA as well as for the entire globe.
Three single-forcing simulations are performed to calculate separately the effects of GHGs, anthropogenic aerosol emissions, and natural forcing on the TAS variations. The effects of interactions between these external forcings on TAS are estimated as the residual term between HST and the sum of the single-forcing runs. Analyses of the simulated data show that the TAS trend in the three periods can be related to distinct effects of the variations in the anthropogenic (GHGs and aerosols) and natural (solar and volcanoes) forcings, as well as their interactions, within the corresponding period. The main results from the analysis of the experiments include: • For P1, the positive TAS trend over the EA land surface is primarily due to the increase GHGs and natural variability. The warming trend by the natural variability is mainly by the long-tem recovery of TAS from the large volcanic activities in late 19th/early 20th century that lasted until 1940. The positive effect of the solar activity is much smaller than the long-term recovery from the volcanic activities. Emissions of anthropogenic aerosols and the interaction term exert little impacts on the TAS trend in the early 20th century period.

•
The negative TAS trend over EA in P2, 1941-1970, results from the combined cooling effect of the increased anthropogenic aerosols and natural forcing as they exceed the combined warming effect of the increased GHG concentration and the interaction term. The anthropogenic aerosol effects in P2 are mainly shaped by the increase in the total aerosol optical depth. The natural forcing effects are due to strong volcanic eruptions during the late P2 period.

•
For P3, the increase in GHGs mainly contributes to the positive TAS trend. The aerosol effects on SW occur mainly in the clear-sky condition. Thus, aerosol direct effects appear to be critical in shaping the aerosol effects on TAS over EA. However, the effects of negative SW of anthropogenic aerosols is overcome by the positive effects of GHG increases to contribute to the positive TAS trend in the late 20th century period. • For P1 and P3, the surface SW forcing is not the driving force behind the TAS change in HST as the results show that the warming effects of the surface LW forcing dominate the cooling effects of the surface SW forcing. Further details of the relationship between the LW forcing and TAS cannot be analyzed because the LW forcing to TAS cannot be separated from the LW response to TAS using the data available to us.

•
It is noticed that the nonlinear interaction effect between the anthropogenic and natural agents substantially contributes to the TAS variations since the mid-20th century (in P2 and P3). The interaction term also contributes to the hydrological variables. The decreased (increased) water vapor and cloud water contents by the interaction term result in the suppression (enhancement) of convective precipitation in P2 (P3).

•
Surface downwelling LW flux is well correlated with the changes in TAS in all experiments and for all time periods. In fact, for P1 and P3, the positive effects of the downwelling LW exceed the negative effects of downwelling SW to contribute to the positive TAS trend in the corresponding period. Quantification of the effects of LW on the TAS trend is not feasible because the simulated LW trends represent the combined effects of the LW forcing and the LW response to the altered TAS.
This study is among the few climate modeling studies that explored the effects of the forcing by natural elements, aerosols, and GHGs in East Asia within a transient framework. The CMIP5 GCM simulations analyzed appears to successfully reproduce the historical TAS in the 20th century and attribute key multi-decadal variability to the historical emissions of GHGs, anthropogenic aerosols, and natural agents, albeit qualitatively. Results in this study are consistent with the previous study [35] that demonstrated that the positive TAS trend in the late 20th century was mostly due to GHG increases, although it was offset significantly by the cooling due to anthropogenic aerosols and the effect of natural forcing was insignificant. However, this study finds that in detail the main forcings on the TAS trends over EA are not uniform but vary according to the analysis periods. In addition, it is also found that the nonlinear interaction among the anthropogenic GHG, aerosols, and natural forcing becomes considerably important in shaping the historical TAS since mid-20th century. Analyses of the effects of the interaction term on the fields related to the atmospheric water cycle suggest that the interaction term affects TAS largely by modifying atmospheric hydrometeors. Detailed mechanisms involved in the interaction term need further research but are beyond the scope of this study. Despite the success, the limitations of this study should be noted. We adopted CMIP5 multi model ensemble analysis to understand the impact of individual forcings on TAS and their uncertainties. However, the interaction term can still be affected by the differences in the specific components of climate forcing and cloud microphysics used among the CMIP5 models. In particular, the uncertainties can be large due to the differences in the physical processes associated with aerosol and cloud interactions among CMIP5 models. Another caveat of this study is the lack of the analysis of the cloud fields, especially the vertical profiles of cloud optical properties. Analyses of the cloud-related fields require special sampling of model fields that is not included in the CMIP5 designs. Uncertainties in attributing the TAS change to forcing terms also arise from our inability to analyze other key forcing terms such as the contribution of the temperature advection related to the altered large-scale circulation because we could not obtain the model data for these analyses. These shortcomings in this study will be explored in future studies.