Multi-model event attribution of the summer 2013 heat wave in Korea Weather and Climate Extremes

To assess the anthropogenic in ﬂ uence on the summer 2013 heat wave in Korea, this study employed a fraction of attributable risk (FAR) approach to three Atmospheric General Circulation Models (AGCMs) with a large ensemble simulation, participating in the C20C þ Detection and Attribution Project. Monthly and daily temper- atures were compared between two experiments. The real world (ALL) experiments were simulated under the observed variations in sea surface temperature, sea ice, greenhouse gas, and aerosol concentrations, while the counterfactual world (NAT) experiments were performed under adjusted boundary conditions by removing anthropogenic warming and with preindustrial levels of greenhouse gases and aerosols. Results from the three AGCMs consistently show that anthropogenic in ﬂ uences had an important role in the extreme heat event over Korea, increasing the chance of the occurrence of extreme warming in summer mean temperature as observed in 2013 by at least 20 times, which supports results from the Coupled Model Intercomparison Project Phase 5 (CMIP5) coupled GCMs (CGCMs). A comparison of individual CMIP5 CGCMs suggests that inter-model difference in FAR values is highly correlated with the amplitude of surface warming centered over Korea, which is also supported by the three AGCMs. Further analysis of individual forcing experiments suggests that the inter-model difference in the regional surface warming is closely linked to the model's response to the aerosol forcing, with stronger in ﬂ uence than that of greenhouse gas forcing. Anthropogenic in ﬂ uences also result in a 5 – 6 times greater likelihood of extreme daily heat events as observed in 2013, which supports a robust mean-extreme relation in the attribution of extreme heat waves. Generally good agreement between AGCM and CGCM results increases the robustness of the conclusion of anthropogenic in ﬂ uences on the summer 2013 Korean heat wave.


Introduction
East Asia (i.e., Korea, eastern China, and Japan) experienced a severe heat wave in summer 2013, which had major impacts on society and the economy. Previous studies compared general circulation model (GCM) outputs simulated with and without human influences and consistently argued that human activities were the main contribution to this extreme event occurrence (Table 1). In South Korea, the summer mean daily minimum temperature in 2013 broke the existing maximum record (Min et al., 2014). To assess the possible impacts of anthropogenic influences on the 2013 Korean heat extremes, Min et al. (2014) compared CMIP5 coupled GCM (CGCMs) with and without anthropogenic forcings using a large-scale sea surface temperature (SST) based indicator that is closely linked with Korean temperature. They found that the summer 2013 like extreme temperature in Korea became 10 times more probable due to human-induced warming.
For eastern China, Zhou et al. (2014) and Ma et al. (2017) found that the anthropogenic influence increased the chance of the 2013 heat wave by 2-3 times based on CMIP5 CGCMs by comparing distributions of anomaly temperatures from simulations performed with and without human influences for a long-term period of 1900and 1955, respectively. Sun et al. (2014 estimated that human activity contributed a greater than 60-times increase in the likelihood of the occurrence of the 2013 hottest summer over eastern China by constructing distributions including future projections of CMIP5 CGCMs. For Japan, anthropogenic influence was assessed to be the main cause with 7-20 times increase in risk of the 2013 heat wave based on MIROC5 AGCM ensemble simulations . Most previous studies of event attribution were based on limited number of CGCMs and the estimated changes in risk of the occurrence of extreme events exhibited a large range from 2 to 60 times, depending on the analysis domains, target variables, climate models, and sampling methods for constructing temperature distributions for the real and counterfactual worlds Min et al., 2014;Sun et al., 2014;Zhou et al., 2014; see Table 1). Moreover, many CMIP5 CGCMs have a relatively coarse horizontal resolution (larger than 100 km), and the robustness of the attribution results needs to be confirmed by comparison with those from other models with higher resolutions. In this respect, the Climate of the 20th Century Plus Project (C20Cþ) Detection and Attribution (D&A) subproject (http://portal.nersc.gov/c20c) has produced a large pool of outputs from atmospheric GCMs (AGCMs) with a relatively high spatial resolution to help understand changes in extreme weather events in the context of past and current climate change (Stone et al., 2018a). Recently, Ma et al. (2017) analyzed large ensembles of two AGCMs, CAM5.1 and MIROC5, participating in this subproject. They suggested that a hot summer over central eastern China such as the 2013 event was 17 times (CAM5.1) and 4 times (MIROC5) more probable due to human influences, confirming the anthropogenic influences but also indicating a strong model dependence of the attribution results.
The C20Cþ D&A project includes the generation and comparison of two climate change scenarios, the observed boundary condition (ALL) and a counterfactual "natural" world (NAT). A large ensemble of single AGCMs has the advantage of generating many simulations to enable robust sampling under different initial conditions. Moreover, model uncertainty can be reduced based on the prescribed observed boundary condition (Christidis and Stott, 2014). It should, however, be noted that a lack of air-sea coupling can affect event attribution results (Dong et al., 2017). Although attribution of large-scale surface air temperature (SAT) changes was suggested to be relatively insensitive to the air-sea coupling (Dong et al., 2017), the impact of air-sea interaction on local scale phenomena such as Korean heat waves needs to be assessed. For this purpose, a comparison of AGCM-based results with those from CGCMs provides an important way of evaluating the robustness of the attribution statements for the local temperature extreme events.
Hence, this study assessed the event attribution of the summer 2013 heat wave in Korea by comparing the C20Cþ D&A models with a CMIP5 multi-model ensemble (MME). The anthropogenic contribution to the extreme event was quantified by comparing the probability of exceeding the observed temperature between simulations with natural forcings alone and simulations with both natural and anthropogenic forcings by using a fraction of attributable risk (FAR) approach. Moreover, we conducted sensitivity tests of the FAR values, focusing on the role of the climate model sensitivity and the difference in the boundary conditions in order to explore physical mechanisms for determining the uncertainty in the event attribution results for the local temperature extremes.

Data and methods
For the observation data, we used daily mean/minimum temperature data measured by 12 Korea Meteorological Administration weather stations from 1954 to 2013. In addition, we used monthly SAT data from HadCRUT4 (Morice et al., 2012) to identify a large-scale indicator associated with Korean summer temperature changes, using an upscaling approach (Min et al., 2014(Min et al., , 2015a. This approach is suitable for analyzing local climate changes, such as those on the Korean Peninsula, because global climate models cannot well capture local-scale weather and climate processes due to their relatively low spatial resolutions. We used three AGCMs that participated in the C20Cþ D&A project: CAM5.1 (Neale et al., 2012;, MIROC5 , and HadAM3P-N96 (Wolski et al., 2014), which were run at resolutions of~1 ,~1.4 , and~1.8 , respectively. The three AGCMs generated two different types of ensemble simulations with and without the effect of anthropogenic influences. The simulations with anthropogenic influences prescribed the observed boundary conditions as greenhouse gases (GHGs), aerosols, SST, sea ice coverage, and land use/cover, referred to as "ALL". The simulation without anthropogenic influences used a CMIP5 estimate of the change in SST, which was subtracted from the observed boundary SST (with sea ice coverage data modified accordingly), which is referred to as "NAT" (Stone and Pall, 2018). GHGs, aerosols, and ozone were set to pre-industrial levels for the NAT simulation (Stone et al., 2018a; the data can be accessed at http:// portal.nersc.gov/c20c/data.html).
The CAM5.1 and HadAM3P-N96 used monthly SST and sea ice coverage from the Hurrell et al. (2008) dataset. The MIROC5 prescribed monthly SST and sea ice coverage which were taken from the HadISST dataset (Rayner et al., 2003). MIROC5 used the prescribed aerosol precursor emissions (i.e., atmospheric chemistry interactions), while CAM5.1 and HadAM3P-N96 used the prescribed aerosol burdens (i.e., black carbon, organic carbon, sulfate, and sea salt) without the chemistry interaction; HadAM3P-N96 used climatological aerosol burdens that were not altered for the NAT simulations. For the event attribution assessment, we used the maximum number of available simulations from each AGCM. The CAM5.1 and MIROC5 runs were trimmed to cover the 2007-2013 period with 100 ALL and NAT ensemble members. The HadAM3P-N96 run was trimmed from 2007 to 2012, and consisted of 50 ALL and NAT members. To assess the ability of the climate models to reproduce the observed long-term interannual variability, we used the ALL simulation for 50 runs from CAM5.1, 10 runs from MIROC5, and 10 runs from HadAM3P-N96 over the 1961-2013 (2012) period. As Chen and Zhou (2017) showed, local and remote SST variability can affect interannual variability of East Asian heat waves. To examine possible influence of SST condition in 2013, we have examined PDF and FAR results using data from the year 2013 for the two AGCMs (CAM5.1 and MIROC5). Results remain unaffected, although slight differences in the probability of occurrence exist (not shown). The influence of SST variability on a given year should be very similar between the ALL and NAT simulations because same SST patterns were prescribed with the only difference in means (delta-SST).
To investigate the dependence of the attribution results on the individual model and/or experimental design, the monthly SATs from the historical, historicalNat, historicalGHG, Representative Concentration Pathway (RCP) 4.5, and preindustrial control (piCTL) experiments from CMIP5 MME (Taylor et al., 2012) were analyzed (Table S1). The ALL simulations were constructed by using data from the historical experiment forced by natural (i.e., changes in solar and volcanic activities) and anthropogenic forcing (mainly changes in GHGs and aerosols) for 1961-2005 and extending them to 2012 using the RCP 4.5 scenario simulations. The NAT simulations for 1961-2012 were prepared from the historicalNAT with natural forcing alone. Similarly, the GHG runs were obtained from historicalGHG experiment integrated with only well-mixed greenhouse gases. We selected 10 CMIP5 models (39 ensembles) that were available for both the ALL and NAT experiments. The CMIP5 models were trimmed from 2007 to 2012 to carry out the event attribution analysis of the summer 2013 extreme event in Korea.
To identify a large-scale SAT indicator of local temperature in Korea, we obtained a correlation map between Korean summer mean daily minimum temperature (Tmin) and East Asian SAT for 1961-2013 using the observations, both of which were June, July, and August (JJA) averages (Fig. 1a). Then, we chose a latitude-longitude box (30-45 N, 115-140 E) with a strong correlation (r > 0.6), hereinafter referred to as NEA (northern East Asia), and calculated the NEA area-averaged SAT anomalies from the observations, three AGCMs, and all CMIP5 simulations (referred to as SAT anomaly). To generate an indicator of summer heat in Korea we used Tmin instead of summer mean daily maximum temperature (Tmax), because Tmin better represents the observed 2013 heat wave with strong intensity and larger spatial coverage than Tmax as reported in Min et al. (2014). Results based on Tmax were found to be similar to Tmin-based results with strong positive correlations with SAT (not shown). The SAT anomalies were calculated relative to the 1971-2000 climatology. For models, climatology was defined as each model's ensemble mean of ALL simulations for 1971-2000 using available runs (see Table S1).
To quantitatively assess the occurrence probability of the extreme event under anthropogenic forcing, we calculated the fraction of  Stone and Allen, 2005). The FAR approach compares the probability of an extreme event occurring between a real world (with human influence) and a counterfactual world (without human influence). The FAR is formulated as FAR ¼ 1 -(P N /P A ), where P N denotes the probability of exceeding the observed event occurring under natural unforced conditions (NAT) and P A represents the same probability estimated under the anthropogenically forced conditions (ALL). The 5-95% uncertainty range of the FAR was estimated from 1000 bootstrap resampling as follows. For each AGCM or CMIP5 MME data, SAT anomalies were first fitted to the kernel distribution (using Gaussian kernel function). Random samples of SAT anomalies were then drawn from the fitted distribution with the same sample size (e.g., 100 for CAM5.1 ALL and NAT) and the corresponding risk ratio (RR ¼ P A /P N ) was calculated. This process was repeated 1000 times, and 5th and 95th percentiles of log (RR) were estimated from 1000 RR values. Here the "basic bootstrap" method (Paciorek et al., 2018) was employed to estimate the confidence interval of log (RR). In this method, the 5th and 95th percentiles are reversed such that confidence interval is estimated as [p m -(p 95 -p m ), p m -(p 5 -p m )] where p 5 , p 95 , and p m are 5th percentile, 95th percentile, and the original all-data estimate of log (RR), respectively. This is to consider possible non-symmetric nature in the tail of the log (RR) distribution although the log (RR) distribution is generally symmetric unlike the FAR distribution (Paciorek et al., 2018). Note that we use this method only when less than 5% of 1000 samples have infinities of log (RR). Finally, the 90% confidence interval of the FAR was obtained by converting the log (RR) percentiles into FAR values.

FAR analysis for the JJA mean temperature
The observed large-scale SAT pattern associated with the summer mean Tmin in Korea was characterized by positive anomalies over the NEA domain, with a maximum over the South Sea and the East Sea/Sea of Japan (Fig. 1a). The time series of the observed NEA averaged SAT anomaly exhibits strong correlation (r ¼ 0.92) with the Korean Tmin (Fig. 1b). The SAT anomaly in 2013 is the second highest after 1994 but with similar amplitude. The three AGCMs and CMIP5 models effectively reproduced the observed large-scale relation between the Korean Tmin and NEA SAT anomalies with correlation coefficients of 0.6-0.98 (not shown). This indicated that the NEA SAT anomaly explains most of the fluctuations in Korean summer temperature in both the observations and model simulations, and can be used as a good indicator of Korean summer temperature. Fig. 1c-f represents the time series of the modeled SAT anomalies averaged over the NEA from the three AGCMs and CMIP5. The ensemble mean of the SAT anomaly was significantly correlated with the observed SAT anomaly, with a correlation coefficient of 0.92 for CAM5.1, 0.90 for MIROC5, and 0.89 for HadAM3P-N96. The interannual variabilities of the SAT anomaly were reproduced well by the three AGCMs, but with some under-estimations when compared with the station observations. The standard deviation (SD) of the NEA averaged detrended SAT anomaly was 0.51 C, which can be larger, partly due to the use of a relatively small number of point measurements, possibly inflating the variance. The ensemble mean of the SD of the detrended SAT anomaly from CAM5.1 was 0.43 C with a 5th to 95th percentile range of 0.37-0.49 C. Similar results were obtained from MIROC5 and HadAM3P-N96, which were 0.42 C (0.34-0.49 C) and 0.44 C (0.40-0.51 C), respectively. In CMIP5, the correlation coefficient (r ¼ 0.49) between the observed SAT anomaly and multi-model ensemble mean (MME) SAT anomaly was less than those of the three AGCMs, representing freely driven SST fluctuations in CGCMs. However, the CMIP5 inter-model spread of the SAT anomaly SD (0.40-0.62 C) is larger than those from AGCM, well covering the observed value, since it includes additional oceanic variability.
To conduct the attribution analysis of the 2013 heat event, we compared the probability of the occurrence of the event under the ALL and NAT simulations. Fig. 2 shows the histograms and kernel curves of the NEA SAT anomalies for the three AGCMs and CMIP5. As described above, the histograms and curves for CAM5.1 and MIROC5 were estimated from 7 years (2007-2013) Â 100 ensemble members but those for HadAM3P-N96 were estimated from 6 years (2007-2012) Â 50 ensemble members. CMIP5 was constructed from 6 years (2007-2012) Â 39 ensemble members (10 models). In CAM5.1 and HadAM3P-N96, the probabilities of the occurrence of SAT anomalies warmer than the observed 2013 event were 4.86% and 3.23% in the ALL simulations, respectively, while the event did not occur under the NAT simulations ( Fig. 2a and c and Table 2). Correspondingly, the contribution of anthropogenic forcing to the observed 2013 event became 100% (FAR ¼ 1, Table 2), representing a negligible influence of natural forcing. For MIROC5, the probabilities of the occurrence of an event hotter than the 2013 observations in the ALL and NAT simulations were 1.47% and 0.08%, respectively ( Fig. 2b and Table 2). The resulting FAR value was 0.95, indicating that anthropogenic influences have increased the risk of heat waves, such as the 2013 event in Korea, by 20 times. The difference in the FAR among three AGCMs can be influenced by the difference in the SAT anomaly over the NEA between the ensemble mean of ALL and NAT realizations (hereinafter referred to as "delta-SAT"). When the delta-SAT increases, the chance of exceeding the observed strength in the NAT simulation would decrease, assuming that the spread of the distribution remains the same. Indeed, MIROC5 had a small delta-SAT (0.29 C) compared with those of CAM5.1 (0.47 C) and HadAM3P-N96 (0.91 C). Considering that prescribed SST was almost  identical for all three models, the SAT differences should be particularly large over land. Indeed, MIROC5 shows a negative delta-SAT over the land area whereas HadAM3P-N96 exhibits larger positive SAT (not shown). Because MIROC5 simulated the aerosol distribution from prescribed aerosol emissions, hot extremes of the ALL simulations from MIROC5 could be more strongly influenced by aerosol cooling feedback Ma et al., 2017). Also, since HadAM3P-N96 retained ALL aerosols in its NAT simulations, the larger positive delta-SAT in HadAM3P-N96 supports the importance of aerosols. Further, Ma et al. (2017) suggested that the difference in the FAR could be explained by the difference in climate sensitivity. The CAM5.1 (2.3 C) and HadAM3P-N96 (2.0 C) models have a higher transient climate response (global and annual mean temperature response at the doubling of CO 2 concentration, estimated from 1% per year CO 2 increase experiment; see below) than MIROC5 (1.5 C). The three AGCMs are, however, not enough to assess the link between FAR and model's climate sensitivity, so we further examined relations between FAR and delta-SAT, climate sensitivity, and the spread of the distribution using CMIP5 CGCMs below (section 3.2). We also conducted the attribution analysis for the 2013 heat wave in Korea using the CMIP5 MME (Fig. 2d) using the same method as done for the AGCMs. The probability experiencing SAT anomalies stronger than the observed 2013 event was 13.83% in the ALL simulation and 2.68% in the NAT simulation ( Table 2). As a result, the FAR value was 0.81, indicative of a 5-times increase in risk due to human influences. Given a delta-SAT of CMIP5 MME (0.66 C) within the range of three AGCMs (0.29-0.91 C), the smaller FAR in CMIP5 MME than in AGCMs seems to be in part due to its larger multi-model spread (see below). The FAR for CMIP5 was a bit smaller than that determined by Min et al. (2014), who found FAR ¼ 0.9, representing a 10-times increase in risk due to anthropogenic influences from the CMIP5 MME using a large-scale SST indicator. However, it is difficult to directly compare the results from Min et al. (2014) to those in this study due to the different CMIP5 model samples, analysis period, and upscaling variable and domain selected (Table 1). When the same CMIP5 model samples were used as in Min et al. (2014), i.e., ALL (31 models -105 ensemble members) and NAT (8 models -38 runs), the FAR was slightly reduced to 0.71. When using the same long-term period to construct the distribution as in Min et al. (2014), the FAR decreased to 0.41 due to the reduction in delta-SAT. Finally, when the upscaling area was extended toward the Fig. 4. Scatter plots between the FAR (y-axis) and sensitivity factors (x-axis): (a) delta-SAT, (b) SD, (c) delta-SAT divided by SD, and (d) TCR. The black dots indicate the individual CMIP5 models, and colored circles indicate three AGCMs (purple: CAM5.1, blue: MIROC5, and orange: HadAM3P-N96) and CMIP5 MME (green). SD values were estimated from NAT for AGCMs and from piCTL for CMIP5 models. The correlation coefficients were calculated using three AGCMs and nine CMIP5 models combined, except for HadGEM2-ES which has negative delta-SAT and FAR. The correlation coefficients in parentheses were calculated using only nine CMIP5 models. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) East Asia domain, a hot summer, such as the 2013 event, could not occur under NAT forcing, and the FAR increased to 1.0. Therefore, the upscaling domain seems to be one of the major factors for the differences in our results with those of Min et al. (2014).

Sensitivity test
We assessed the sensitivity of FAR for the summer 2013 heat wave event in Korea focusing on the roles of delta-SAT, the spread of distribution, and the models' climate sensitivity. For the individual CMIP5 model, relatively small ensemble size for ALL and NAT runs (Table S1) makes it difficult to calculate the FAR for the individual model. Therefore, assuming that the CMIP5 individual models and CMIP5 MME have the same shape of the SAT anomaly distributions, we calculated the FAR of the CMIP5 individual model after adjusting the CMIP5 MME distribution to the mean and spread of the individual models. First, for each model, means are obtained from the ALL and NAT simulations during the 2007-2012 period and SD is estimated from the piCTL experiment after removing climate drift for the available whole period. Next, CMIP5 MME distribution is normalized for ALL and NAT simulations respectively, by removing its mean and dividing it with its SD. Finally, we multiply the normalized CMIP5 distributions by the SD of the individual model and add the corresponding means for ALL and NAT from the model. Fig. 3 shows the histogram and kernel density distributions for the four models obtained by adjusting the CMIP5 MME distribution as described above. The difference in the mean SAT between ALL (0.58 C) and NAT (À0.08 C) from CMIP5 MME was 0.66 C, and the SDs during 2007-2012 were 0.53 C (ALL) and 0.59 C (NAT), respectively. For the FAR sensitivity test, we used one SD (0.51 C, vertical dashed line in Fig. 3) of the observed SAT anomaly as a threshold for FAR calculations because, in several CMIP5 models, extremes corresponding to the 2013 SAT anomaly (1.20 C, vertical solid line in Fig. 3) did not occur in NAT runs, making it hard to assess inter-model FAR differences (Fig. S1). The FAR of the CMIP5 MME for one SD of the observed SAT anomaly was 0.71. As an example, for CanESM2, the delta-SAT (1.24 C) was larger than that of the CMIP5 MME and the spread (SD ¼ 0.37 C) estimated from piCTL was smaller than that of the CMIP5 MME. After adjusting CMIP5 MME with the CanESM2 delta-SAT and spread, the FAR increased to 0.96. In contrast, the CSIRO-Mk3-6-0 model had a smaller delta-SAT than that of the CMIP5 MME and a similar spread to the CMIP5 MME, resulting in a decrease in FAR (0.37). When using GISS-E2-R with a similar delta-SAT to CMIP5 MME, the FAR of the adjusted CMIP5 MME was 0.82. The IPSL-CM5A-LR had the largest delta-SAT (1.71 C), and the SD was similar to that of the CMIP5 MME; the FAR was estimated as 0.99. These examples suggest important role of delta-SAT in determining FAR.
To further check the delta-SAT and FAR relationship, a scatter plot is drawn using the CMIP5 individual models which were obtained from adjusting CMIP5 MME distributions as explained above (Fig. 4a, black  dots), and the three AGCMs (color dots). Green dots represent CMIP5 MME. Delta-SAT values (ALL minus NAT runs) are all positive. One exception is HadGEM2-ES model which has a colder surface condition in ALL than in NAT, which is likely due to too strong response to the aerosol forcing (see below). The FAR values of the three AGCMs (color dots) fell on the CMIP5 relationship line, indicating that the differences in the FAR from three AGCMs can be explained by the delta-SAT. Indeed, the FAR and delta-SAT values were significantly correlated across all CGCM and AGCM models (r ¼ 0.86) at the 5% significance level, which remains very similar when using CGCM models only (r ¼ 0.85).
We examined the influence of the spread of the distribution, which was defined by the SD estimated from the large-ensemble NAT simulations during 2007-2012 (2013) for individual AGCMs. The SD was estimated from piCTL runs for individual CMIP5 models. The linear Fig. 5. Scatter plot of (a) delta-SAT and SAT response to AER, (b) delta-SAT and SAT response to GHG, (c) SAT response to AER and FAR, and (d) SAT response to GHG and FAR using ten CMIP5 models. The inter-model correlation coefficients are calculated, except for HadGEM2-ES. The correlation coefficients in brackets are based on the ten CMIP5 models, including HadGEM2-ES. relationship between FAR and the spread of the distribution was found insignificant, indicating that internal variability did not affect FAR much across models (Fig. 4b). However, this is not inconsistent with Bellprat and Doblas-Reyes (2016) who suggested the importance of the spread in FAR estimation. The influence of the spread on FAR is difficult to identify in our multi-model setting where delta-SAT also varies across models.
To further analyze the relative contributions of delta-SAT (signal) and the spread of the distribution (noise) to FAR, we determined the relationship between the signal-to-noise ratio (SNR) and FAR (Fig. 4c). The FAR and SNR showed a strong positive correlation (r ¼ 0.88) similar to that for the delta-SAT, suggesting the dominant influence of delta-SAT on the attribution statements for Korean heat wave. The FAR of the model, in which the signal was relatively larger than noise, became near 1 despite the low observed strength (i.e., one SD). Also, difference between CMIP5 MME and AGCMs was reduced in the SNR-FAR relation scatter, indicating that the larger inter-model noise in CMIP5 MME indeed made its FAR smaller as discussed above (Fig. 2).
Further, we evaluated the impact of TCR to the FAR, as discussed in Ma et al. (2017). TCR is a useful metric for assessing the climate sensitivity of a GCM, which is defined as the global and annual mean temperature change at the timing of CO 2 doubling following a linear increase in CO 2 concentration (1% per year) over a period of 70 years . The TCR for the individual CMIP5 model, CAM5.1 and MIROC5 was obtained from Flato et al. (2013;Table 9.5) and the TCR for HadAM3P-N96 was assumed to be same as that for HadCM3 (Randall et al., 2007, Table 8.2). Fig. 4d shows the relationship between the TCR and FAR for the CMIP5 individual model and three AGCMs. The relationship is generally positive, meaning that models with higher TCR tend to have larger FAR, which can be seen also from three AGCMs. However, the FAR-TCR relationship is not statistically significant.
Overall, the analysis results of inter-model relation suggest that delta-SAT had a more significant role in determining the FAR for the summer heat wave event in Korea than model's TCR to the greenhouse gas forcing or the spread of distribution. We also tested the sensitivity to the different event thresholds such as 0.5 SD, 1.5 SD, and 2 SD. For example, 2SD (1.02 C) is very close to the previous event occurred in 2010 (1.13 C, see Fig. 1b). Results were very similar to Fig. 4 based on 1 SD threshold (not shown), indicating the robust relationship between FAR and delta-SAT, delta-SAT/SD and TCR.
The next question would be then what determines delta-SAT in this region. For simplicity, we assume that delta-SAT can be approximated to be a combined response to GHG forcing and aerosol forcing. To investigate which forcing component affected delta-SAT more importantly, inter-model relationship was examined between delta-SAT and SAT response to GHG and aerosol forcings using CMIP5 models. Regional SAT response to GHG was estimated from GHG results during 2007-2012 while regional SAT responses to aerosol forcing (AER) was estimated from a residual by subtracting the responses to GHG and NAT results from ALL results during the same period. We note that this estimation of AER as a proxy may not be representative of actual aerosol influence on temperature due to possible influences of other forcing factors such as land use change and ozone forcing (Forster et al., 2013;Shindell et al., 2015).
Results show that GHG and AER exhibited a consistent warming and cooling effect, respectively (Fig. 5). Interestingly, delta-SAT were highly correlated with the mean SAT anomaly response to the aerosol forcing, with correlation coefficients of 0.90. This indicates that models having stronger aerosol cooling responses near South Korea have cooler surface temperature responses in ALL, and hence less difference from NAT results. A strong positive correlation was also found between FAR and aerosol response (r ¼ 0.81, Fig. 5c). Further, we examined the relationship between delta-SAT and SAT response to aerosol forcing using anthropogenic aerosol (AA) only forcing runs available from five CMIP5 models (CanESM2, CSIRO-Mk3-6-0, GISS-E2-H, GISS-E2-R, NorESM1-M) for 2007-2012. Results show positive inter-model correlations between delta-SAT and SAT response to AA forcing and between FAR and SAT response to AA forcing (Fig. S2). Although statistically insignificant, this supports the important role of aerosol forcing in determining FAR spread among models. The weaker correlation from AA runs seems to be due to a small number of models and also possible influences of other anthropogenic forcings on local SAT as discussed above, which warrants further investigation.
Overall, our results suggest that inter-model difference in delta-SAT and FAR over this region is largely determined by model's response to aerosol forcing. This is generally consistent with the previous findings for the global mean and Northern Hemispheric surface temperature trends (Shindell et al., 2015). In contrast, delta-SAT and FAR was not sensitive to GHG response ( Fig. 5b and d), indicating a weak link between the spread in model-derived FAR estimates and the models' TCRs. Although AER contribution to delta-SAT in the three AGCMs is not available, MIROC5 model seems to have a stronger response to the aerosol forcing, as Ma et al. (2017) discussed. The role of aerosol forcing in determining inter-model differences at regional to local scale climate changes requires further investigation using relevant model experiments (e.g., Detection and Attribution Model Intercomparison Project/CMIP6, Gillett et al., 2016).

FAR analysis for the JJA maximum daily temperature
In terms of impact of extremes, seasonal maxima of daily temperature can be more relevant than their seasonal means due to its stronger intensity. During the 2013 summer, Korea experienced a record-breaking summer maximum of daily minimum temperature (TNx). Fig. 6 shows the time series of TNx based on 12 station observations and the three AGCMs. The TNx for the climate model was calculated as the average for the TNx over the Korea area (34-38 N, 125-130 E, see small box in Fig. 1a). The correlation coefficient between the Korean TNx and TNx simulated from CAM5.1 was 0.72, indicating that the model effectively explained observed fluctuations in the summer extreme temperature in Korea. For example, the TNx simulated from CAM5.1 showed a record high in 2013. The amplitude of the interannual variability of TNx was also reproduced well by CAM5.1. The ensemble mean of the SDs of the detrended TNx simulated from CAM5.1 was 0.62 with a 5th to 95th percentile range of 0.54-0.74, covering the observed value (0.7 C). MIROC5 reasonably simulated the interannual variability in TNx, with a correlation coefficient of 0.60, whereas the ensemble spread of TNx was smaller with 90% range of 0.45-0.58 C. HadAM3P-N96 had a larger fluctuation of TNx interannual variability, with a 5th to 95th percentile range of 0.74-1.03, due to a positive-skewed distribution (see below).
The distributions of TNx simulated from the three AGCMs of the ALL (2007)(2008)(2009)(2010)(2011)(2012) and NAT (2007NAT ( -2012 runs were compared with the observed 2013 TNx (Fig. 7). In CAM5.1, events with TNx values stronger than the observed 2013 value were extremely rare in the NAT runs (0.58%). The chance of a 2013-like extreme event increased to 2.86% in the ALL runs. The corresponding FAR value was 0.80 (90% range: 0.27-0.91), indicating a 5-times increase in likelihood due to human influences (Table 2). In MIROC5, the probability of exceeding the observed 2013 TNx from the ALL and NAT runs was 1.29% and 0%, respectively, indicating that extreme temperatures such as those observed in the 2013 event have never occurred in the absence of anthropogenic climate change. The probability simulated from HadAM3P-N96 for years more extreme than 2013 in the ALL and NAT runs was 9.59% and 1.77% (Fig. 7c), corresponding to a FAR of 0.82 (90% range: 0.47-0.90). These results based on daily temperatures indicate that anthropogenic influences have increased the risk of extreme daily heat events, such as that in 2013, in Korea by 5-6 times (Table 2), generally consistent with the results for the summer mean temperature. As shown in the scatter plot in Fig. 7, the summer mean temperature (xaxis) was correlated with extreme daily minimum temperature (y-axis) in both the ALL and NAT simulations, with ANT distributions being shifted from NAT distributions to the upper right direction, indicating warming of both seasonal mean and seasonal extreme temperatures during summer.
When computing the FAR corresponding to one SD (0.7 C) for the observed TNx, the FAR for three AGCMs was 0.68 (CAM5.1), 0.76 (MIROC5), and 0.79 (HadAM3P-N96), respectively. We found above that the delta-SAT had an important role in determining the FAR for JJA mean temperature based on the nine CMIP5 models (Fig. 4). Based on the three AGCMs, the FAR values for the summer maxima of daily minimum temperature also seem to be related to delta-SAT such that HadAM3P-N96 with the larger delta-SAT (1.17 C) has larger FAR. Considering that response of daily temperatures will be more uncertain, more affected by other local factors like synoptic conditions, further investigation of the associated physical mechanisms is needed in the future work (Min et al. 2015b).

FAR sensitivity to the threshold
It is useful to check FAR values to different temperature thresholds that can occur in the future (e.g., Christidis et al., 2015). For this purpose, here we calculated FAR values for hypothetical temperature anomalies from À3 C to þ3 C in 0.1 C intervals for the JJA mean temperature and TNx. Fig. 8 displays the resulting FAR curves for each AGCM. The slopes of the FAR curve for the three AGCMs were similar, but HadAM3P-N96 showed the fastest saturation of FAR values to 1 at SAT anomalies around þ0.5 C (Fig. 8a). The FAR values of CAM5.1 and MIROC5 were saturated to 1 at higher SAT anomalies of þ1.0 C and þ1.3 C, respectively. The faster increase of FAR in HadAM3P-N96 seems to be because its NAT simulations retain anthropogenic levels of aerosols. The saturation time seems to be related to delta-SAT. The delta-SAT for HadAM3P-N96 was the largest among the three AGCMs (Fig. 2c). The FAR curves of the three AGCMs became saturated to 1 by the time JJA mean SAT anomalies increased to 1.3 C (Fig. 8a), indicating that any SAT anomaly larger than 1.3 C could not occur without anthropogenic influences. The FAR curve of the CMIP5 increased at a slower rate than those of the three AGCMs, which is likely due to the wider spread of the CMIP5 MME distributions by including additional endogenous variabilities in ocean temperatures, as discussed above (Fig. 2).
For TNx, the occurrence time of the risk was the fastest in HadAM3P-N96, similar to the JJA mean SAT anomaly (Fig. 8b). The FAR curve of MIROC5 reached 1 fastest, whereas that of HadAM3P-N96 reached 1 last. This trend seems to be partly related to the distribution shape. The distributions of TNx for both the ALL and NAT runs from HadAM3P-N96 exhibited positive skews with longer tails (Fig. 7c), which is consistent with the larger year-to-year variability (Fig. 6c). In contrast, the distributions of TNx from MIROC5 showed smaller differences in the mean and spread than the other AGCMs, which induces a stronger slope of FAR increase and the faster saturation. CAM5.1 is in between the other two models.

Summary and conclusions
This study conducted the event attribution analysis of the 2013 summer heat wave in Korea by using a large ensemble of three AGCMs (CAM5.1, MIROC5, and HadAM3P-N96) that were included in the C20Cþ D&A project in comparison with CMIP5 MME. Based on the FAR approach, the probability of the occurrence of an extreme event was compared between simulations with natural forcings alone and simulations with both anthropogenic and natural forcings. The comparison revealed that the unusual heat wave in Korea could not be explained without the inclusion of anthropogenic influences, with a 20-times increase in likelihood estimated from MIROC5. Anthropogenic influences were found to account for 100% as estimated from CAM5.1 and HadAM3P-N96. These largely supported the CMIP5 MME based results showing 5-times increase in probability. Similar results between AGCMs and CGCMs further indicate that the attribution conclusion for local-scale heat waves like the 2013 hot summer in Korea is not affected much by air-sea interaction, supporting Dong et al. (2017) who suggested relative insensitiveness of the large-scale SAT attribution to the air-sea coupling.
When examining the CMIP5 individual models, the inter-model uncertainty in the FAR for the heat event was found to be more sensitive to inter-model difference in the delta-SAT than those in the climate sensitivity to GHG forcing and the spread of distribution. The three AGCMs also followed the CMIP5 inter-model relation between FAR and delta-SAT. We further explored possible link of aerosol forcing with the inter-model spread in delta-SAT over Korea. Results suggested that intermodel spread in delta-SAT (and hence FAR) was largely explained by the model's response to the aerosol forcing, with more importance than the model's response to the greenhouse gas forcing over this region. It should be noted that the estimation of the cooling response to the aerosol forcing was computed as a residual of non-GHG anthropogenic forcing from available model simulations (ALL -NAT -GHG). Results using anthropogenic aerosol forcing only runs from limited number of models largely support our conclusions, but the high sensitivity of aerosol influences on FAR warrants further investigation. Using the three AGCMs, this study also performed an attribution assessment for the summer maximum daily minimum temperature (TNx) over South Korea, which broke the record for warmest value in Korea in 2013. The results showed that anthropogenic influences resulted in a 5-6-times more likely occurrence of daily extreme temperature, consistent with the results based on the JJA mean temperature. In conclusion, the results of this study support previous findings that the 2013 extreme heat event in Korea is attributable to the anthropogenic forcing, but the model's response to aerosol forcing may impact the significance of the attribution results.