Advanced risk-based event attribution for heavy regional rainfall events

Risk-based event attribution (EA) science involves probabilistically estimating alterations of the likelihoods of particular weather events, such as heat waves and heavy rainfall, owing to global warming, and has been considered as an effective approach with regard to climate change adaptation. However, risk-based EA for heavy rain events remains challenging because, unlike extreme temperature events, which often have a scale of thousands of kilometres, heavy rainfall occurrences depend on mesoscale rainfall systems and regional geographies that cannot be resolved using general circulation models (GCMs) that are currently employed for risk-based EA. Herein, we use GCM large-ensemble simulations and high-resolution downscaled products with a 20-km non-hydrostatic regional climate model (RCM), whose boundary conditions are obtained from all available GCM ensemble simulations, to show that anthropogenic warming increased the risk of two record-breaking regional heavy rainfall events in 2017 and 2018 over western Japan. The events are examined from the perspective of rainfall statistics simulated by the RCM and from the perspective of background large-scale circulation fields simulated by the GCM. In the 2017 case, precipitous terrain and a static pressure pattern in the synoptic field helped reduce uncertainty in the dynamical components, whereas in the 2018 case, a static pressure pattern in the synoptic field provided favourable conditions for event occurrence through a moisture increase under warmer climate. These findings show that successful risk-based EA for regional extreme rainfall relies on the degree to which uncertainty induced by the dynamic components is reduced by background conditioning.


INTRODUCTION
It is widely acknowledged that global warming significantly affects the features of extreme events 1 . Several studies have focused on identifying changes in the occurrence of extreme events that could be linked to anthropogenic activities or natural variability. However, quantifying changes in the likelihood of rare events is not straightforward. EA based on historical/counterfactual-nonwarming simulations by state-of-the-art climate models is a possible solution for quantitatively estimating the impact of human activities on a specific event 2 . With regard to EA, two questions have been raised: (1) how has anthropogenic climate change altered the probability of occurrence of individual extreme events (known as the 'risk-based' approach; hereafter risk-based EA) 3,4 and (2) how has the severity (i.e. magnitude or intensity) of extreme events changed due to climate change (known as the 'storyline' approach; hereafter storyline EA) 5 . For the risk-based approach, large-ensemble simulations targeted at one single event are required to represent the probability density function (PDF) of the event occurrence. However, for the storyline approach, model boundary conditions are strictly constrained by realistic background atmospheric conditions to ensure the occurrence of the event of interest. Through these approaches, a reasonable degree of consensus can be obtained on the attribution of temperature extremes to anthropogenic climate change, e.g. refs. 6,7 . However, there is less agreement in risk-based EA regarding extreme precipitation [8][9][10][11][12][13][14][15] . Globally, the moisture increase caused by warmer atmospheric temperature agrees with the Clausius-Clapeyron relation of~7% saturated vapour per Kelvin 1 ; however, local precipitation changes are influenced by moisture availability, atmospheric circulation and vertical stability 16,17 . Furthermore, disastrous extreme precipitation is often strongly linked to mesoscale rainfall systems and specific geography 18 . Precipitation takes different forms of atmospheric circulation in each location. The conventional global climate models used for risk-based EA typically cannot capture these features of extreme precipitation events 9 ; therefore, a high-resolution regional model is required. Several previous risk-based EA studies incorporated regional climate models 11,12,14,15 , although their model evaluations were primarily based on heavy rainfall statistics, and representation of background synoptic circulation patterns and mesoscale rainfall systems were not fully considered.
Thus far, high-resolution RCMs have been mainly used in storyline EA. Kawase et al. 19 simulated the heavy rain event of July 2018 in Japan using an RCM with realistic initial and boundary conditions based on the reanalysis data and found that the record-breaking precipitation could be attributed to recent rapid warming through comparison between hindcast and detrended experiments. This technique is similar to the 'pseudo-global warming' approach and has been extensively used for storyline EA of tropical cyclones (TCs). Previous studies conducted a storyline EA aimed at TC-induced extreme events based on RCM and variable-resolution GCM simulations and showed significant impacts of global warming on rainfall amounts, wind speed, storm size and storm surge height [20][21][22][23][24] . However, storyline attribution cannot make a statement about the changes in cyclogenesis rates or the probability because it is conditional on the event's underlying existence 20,24 .
Thus, complete assessments of the effect of anthropogenic climate change on extreme rainfall events must be based on a combination of the storyline EA and risk-based EA. Wang et al. 25 1 Meteorological Research Institute, Japan Meteorological Agency, 1-1 Nagamine, Tsukuba, Ibaraki 305-0052, Japan. 2 Atmosphere and Ocean Research Institute, the University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8568, Japan. 3 Japan Agency for Marine-Earth Science and Technology, 3173-25 Showa-machi, Kanazawa-ku, Yokohama, Kanagawa 236-0001, Japan. 4 National Institute for Environmental Studies, 16-2 Onogawa, Tsukuba, Ibaraki 305-8506, Japan. ✉ email: yimada@mri-jma.go.jp incorporated both of these approaches to examine the effects of climate on Hurricane Harvey's extreme rainfall. Their storyline approach was similar to existing storyline EA based on RCMs, whereas their risk-based discussion on event frequency was based on large-ensemble simulations of relatively long-term future projections by a 100-km coupled GCM, where the poor resolution prevented accurate representation of event-specific convective systems. Based on this, large-ensemble simulations for risk-based EA should be conducted using a high-resolution RCM, although huge computational cost is a serious problem.
In this study, we conducted high-resolution large-ensemble simulations using RCM that were downscaled from large-ensemble GCM simulations to achieve risk-based EA of local heavy rainfall. Because extreme rainfall takes on different forms of atmospheric circulation in each location, each event was examined from the viewpoint of rainfall statistics simulated by the RCM and background large-scale circulation fields simulated by the GCM responsible for rainfall systems that differ according to each location. Global models specialise in capturing the background pathway of large-scale circulation changes associated with global warming. Hence, a large-ensemble set of GCM and RCM simulations is necessary to provide information on changing climate extremes with a reasonable understanding of background processes (Fig. 1). Here, we prepared a 100-member ensemble set of historical and non-warming simulations (hereafter HIST and NonW, respectively) with 100 pairs of Meteorological Research Institute atmospheric general circulation model 3.2 (MRI-AGCM3.2) with 60-km-mesh and non-hydrostatic regional climate model (NHRCM) with a 20-kmmesh for East Asia from 1981 to present, for a total of more than 3000 years for each ensemble (see 'Methods') 26 . The 20-km-mesh is not necessarily enough to represent fine rainfall systems, such as meso-γ scale line-shaped precipitation systems, but is sufficient for representing meso-α or meso-β scale systems such as fronts, tropical depressions and squall lines 27 . The difference between the 60-km and 20-km grid spacing is critical to the representation of Japanese weather. The 60-km resolution cannot distinguish the four main Japanese islands, whereas the 20-km resolution can resolve those islands and can distinguish climate features in the windward and leeward side ( Supplementary Fig. 1).
Based on this methodology, three episodes of extreme precipitation in Japan were attributed to anthropogenic climate change. The July 2018 episode (Case2018) was induced by abnormal moisture inflow toward a stationary rainband in the coastal regions of Japan's Inland Sea (Fig. 2a), and was characterised by its relatively long duration and extensive area of impact 28 . The heaviest 72-h accumulated rainfall was recorded at 123 observation sites from 28 June to 8 July, 2018 (reported by the Japan Meteorological Agency, JMA). The July 2017 episode (Case2017) was associated with abnormal moisture inflow toward the west side of the mountain range on the main island of Kyushu (Fig. 2b). The daily maximum rainfall exceeded 500 mm per day in the worst-hit area from 5 to 6 July (JMA). Both cases accompanied the active Baiu rainband. The ascending current and instability of the Baiu rainband were enhanced by the upper-level westerly jet and travelling synoptic waves, the mid-level advection of warm and moist air influenced by the South Asian thermal low, and low-level southerly moisture transport associated with an enhanced north Pacific subtropical high (NPSH) 29 . All of these circulation features were observed in Case2018 30 . For Case2017, enhanced NPSH was notably observed. The July 1993 episode (Case1993) was caused by a series of TC landfalls on the east side of the Kyushu Island mountain range (Fig. 2c). Two TCs, following similar tracks, made landfall 2 days apart from each other. Daily maximum rainfall exceeded 300 mm per day at several observation sites (JMA).

Validation of model reproducibility
First, we confirmed whether these rainfall systems and responsible atmospheric circulations in each location are successfully classified Extremely strong southerly moisture flowed into western Japan due to the zonal dipole of western cyclonic and eastern anticyclonic circulation patterns. The dipole pattern is well captured by the 60-km AGCM simulation. The lower two panels show the accumulated rainfall amount over the same period [mm], where the right was obtained from the JMA meso-reanalysis data (1-km horizontal resolution), and the left is the result of one of the 100-member NHRCM historical simulations. The 20-km NHRCM roughly captured the regional difference in rainfall amount depending on mountain ranges. The middle panel shows the accumulated rainfall output from the 60-km AGCM. The rainy area was diagnosed in accordance with the synoptic convergence zones and regionality over the lands could not be distinguished. by our simulations. This was done by calculating regression coefficients between each local heavy precipitation frequency index calculated from NHRCM and large-scale circulations (850 hPa geopotential height, column moisture flux and TC density) from MRI-AGCM3.2 using 3000-year samples (100 ensemble members from 1981 to 2010). Figure 2e indicates that the heavy rainfall frequency on the west side of the Kyushu mountain range increases when the NPSH is active and southwesterly moisture flow around the NPSH converges in this region. The negative height anomalies around 40°N correspond to the enhanced Baiu rainband. In the upper troposphere, the features of synoptic waves travelling along the westerly jet are also significant ( Supplementary Fig. 2). This indicates that the specific features of the enhanced Baiu rainband identified by Sampe and Xie 29 are well-represented in our model simulations. TCs do not significantly affect heavy rainfall in this region, likely due to the active NPSH (Fig. 2h). The southwesterly moisture inflow is also present over eastern Kyushu, but it rarely causes heavy rainfall on the eastern side, likely due to the Kyushu mountain range being in a northsouth configuration which would blocks the moisture inflow. Instead, the major cause of the heavy rainfall is TCs which approach eastern Kyushu (Fig. 2i), suggesting that NPSH is not the key factor (Fig. 2f). This is also marginally observable when the same regression map was drawn using the site-observation and reanalysis data with less significance due to the smaller sample size ( Supplementary Fig. 3), indicating the advantage of using large-ensemble data. In the relatively dry regions around Japan's Inland Sea (Fig. 2d, g), the situation is similar to TC-induced heavy rainfall in Fig. 2f, i, and the dipole of the western trough and the eastern ridge partly promotes southerly moisture flow, although both relationships are smaller (but statistically significant) than the other two types. This implies that the patterns in Fig. 2d, g are the result of a combination of factors.
We also compared the leading modes of interannual variability in the circulation patterns over East Asia between the MRI-AGCM3.2 ensemble simulations of HIST and the JRA55 reanalysis from 1981 to 2010, confirming that the leading Rossby wave patterns over East Asia are well-represented in the MRI-AGCM3.2 large-ensemble simulations (Supplementary Fig. 4). In addition to these meso-α to meso-β scale rainfall systems, smaller-scale phenomena such as line-shaped precipitation systems and orographic precipitation induced by complex topography are often observed during heavy rainfall events in these regions. As mentioned above, the 20-km NHRCM is insufficient to resolve these meso-γ scale phenomena. Thus, in this study, we considered rainfall systems larger than meso-β scale, which were the key contributors to triggering the extreme rainfall events of focus.
Attribution of long-term trends in heavy rain frequency Figure 3 shows the difference in heavy rainfall frequency (the number of days with more than 100 mm per day) during the month of July from 1981 to 2010 between the HIST and NonW simulations. It was noted that the impact of global warming on the July 1981-2010 heavy rainfall frequency was strong on the western side of the mountain range because the moisture increase due to warming have directly enhanced the southwesterly moisture convergence 18 . However, the difference was relatively small in eastern Kyushu and the regions of Japan's Inland Sea. In the eastern Kyushu region, large uncertainty in simulated TCs prevented accurate detection of significant differences 18 . In the coastal regions of Japan's Inland Sea, heavy rainfall differences were difficult to attribute to climate warming due to various kind of atmospheric variabilities from year to year. Figure 3 implies that Case2018 (Japan's Inland Sea; various factors) and Case1993 (eastern Kyushu; TC-induced) extreme rainfall events are less likely to be attributed to anthropogenic climate change than Case2017 (western Kyushu; orographic moisture convergence). Figure 4 compares Case2018-like, Case2017-like and Case1993-like event probabilities between the HIST and NonW simulations. Three-day accumulated rainfall was used for Case2018, and daily rainfall was used for Case2017 and Case1993. Because the simulated daily rainfall of 20 km-RCM slightly underestimated the frequency of intense rainfall of more than 200 mm per day 18 , we used a 50-year return value of each 3000 sample (100 members from 1981 to 2010) of the HIST simulation as a threshold instead of observed values (see 'Methods'). Contrary to Fig. 3, Fig.  4 shows a significant increase in extreme rainfall frequency not only in Case2017-like events (western Kyushu) but also in Case2018-like events (Japan's Inland Sea) when we focused on each episode. Figure 4a  which are the 2th-98th percentiles of 5000 bootstrap random samples) to 4.82% (4.38-5.27%) due to human activity. In the same way, Fig. 4b (Case2017-like events) indicates that the return period (probability) of daily maximum rainfall changed from 53.5 years (1.87% with a confidence interval of 1.61-2.12%) to 36.0 years (2.78% with a confidence of 2.48-3.08%) due to global warming. Such an influence was not detectable in the TC-induced type events, as in Case1993 (Fig. 4c). We analysed the background specific synoptic pressure patterns in Case2017 and Case2018. The observed and ensemble-mean circulation anomalies are shown in Fig. 5. Anomaly patterns from the reanalysis data and the ensemblemean simulation were partly dissimilar to each other because the anomaly by the JRA55 reanalysis (Fig. 5a, c) was only one sample including stochastic intrinsic variability, whereas the simulated anomaly (Fig. 5b, d) was the average of 100 members where stochastic intrinsic components were cancelled out. The signals obtained from the ensemble-mean values indicate that the background conditioning behind the event of focus is forced by boundary conditions prescribing MRI-AGCM3.2, such as sea surface temperature (SST). The anomalies for Case2017 (western Kyushu) (Fig. 5c, d) show patterns conditioned by the enhanced NPSH and southwesterly moisture flow, which indicates high heavy rainfall potential as shown by Sampe and Xie 29 . Similarly, the observed and ensemble-mean anomalies for Case2018 (Japan's Inland Sea) (Fig. 4a, b) show that the background field was conditioned by the zonal pressure dipole and the southerly moisture flow. This means that, if the 2018 episode was separated from the other events in the coastal region of Japan's Inland Sea, the conditions of Case2018 and Case2017 would be similar in that moisture convergence was the main driver of high heavy rainfall potential. This is why Case2018-like events were also attributable to anthropogenic global warming.

DISCUSSION
In general, EA results largely depend on background conditioning. Local rainfall changes depend on both background circulation (dynamic) and moisture (thermodynamic) changes. Increased water vapour is a common trend due to global warming, but the reasons for changes in dynamic components are still poorly understood and are thought to be region-dependent 31 . We used the pressure dipole indices to characterise the circulation features for Case2017 (see 'Methods'). Regarding background circulation changes, there is little difference between the HIST and NonW simulations (Fig. 4d), although the 2017 PDFs (both from the HIST (red solid line) and NonW (blue solid line) simulations) are completely separated from the PDF of long-term statistics (shadings). Generally, detection of circulation changes is difficult in the current warming stage because of the large natural variability within dynamical processes. Therefore, the existing EA approach experienced difficulty in detecting extreme rainfall changes. The significant findings of this study are that the ensemble simulations represented the 2017 orographic moisture convergence and 2018 stationary moisture advection well in the ensemble-mean fields as signals (Fig. 5b, d), and that the event probability (Fig. 4a, b) were conditioned on the dynamic situation leading to the event. Because of its moderate constraint on the moisture advection, the role of moisture increase due to global warming was noted in both Case2017 and Case2018. The TCinduced extreme rainfall (Case1993), however, was not attributable to anthropogenic climate change because the occurrence of heavy rainfall depended more on the TC approach rather than on the moisture increase due to global warming. Historical changes in TC density are not detectable at this stage because of large natural variability and problems in reproducibility of TCs by the GCM outside the boundary of the RCM. Furthermore, recent studies have shown that aerosol cooling reduces TC potential intensity

Change in heavy rainfall frequency (HIST -NonW, 1981-2010)
days per month Fig. 3 Past changes in heavy rainfall frequency. Changes in heavy rainfall frequency (the number of days with more than 100 mm per day) in July between the HIST and NonW simulations from 1981 to 2010. The heavy rainfall frequency was obtained from the RCM outputs.
Y. Imada et al. and counteracts an increasing effect due to greenhouse gas warming, which can be a cause of the uncertainty of a TC simulation 32,33 .
Our particular interest is that a comparison of the long-term probability (e.g. 1981-2010 as in Fig. 3) between the HIST and NonW ensembles reveals an almost undetectable difference in the coastal region of Japan's Inland Sea. This is thought to be due to the long-term runs, which include various types of rainfall systems with large uncertainties due to natural variability. However, by the model background conditioning of the 2018-like specific pressure patterns, the dynamic circulation uncertainty decreased, and we detected a change in event probability and attribute it to anthropogenic global warming.
We classified the regional rainfall processes according to the large-scale circulation types and connected each process change to anthropogenic global warming by analysing regional model outputs concurrently with parent global model outputs. This study is an attempt to extend the large-ensemble pairs of regional and global model simulations for the purposes of risk-based EA. If the large-ensemble outputs represent dynamic background conditions leading to a target event in their ensemble-mean field, the risk of heavy rainfall events can be confidently attributed to human activity 5 . Regarding Case2018, we also conducted a conventional storyline EA using NHRCM with a 5-km horizontal resolution 19 . By combining the storyline EA with the risk-based EA, we made more comprehensive assessments of the effect of anthropogenic climate change on extreme rainfall events. This study offers a new direction in the field of high-resolution RCM studies.

The observational dataset
The station dataset for daily precipitation was obtained from the Automated Meteorological Data Acquisition System (AMeDAS) maintained by the Japan Meteorological Agency. Approximately 1300 stations at~17km intervals were available. The analysis period from 1981 to the present day used in this study was determined by the availability of AMeDAS stations, which began surface observations around 1980. The observed anomaly was defined with a reference period between 1981 and 2010. For geopotential heights and column-accumulated water vapour fluxes, we used the 55-year Japanese Reanalysis (JRA55) dataset (ref. 34 ). Typhoon best track data were obtained from the Regional Specialized Meteorological Centres Tokyo -Typhoon Centre.

Large-ensemble simulations
Our large-ensemble pairs of GCM and RCM simulations were part of the Database for Policy Decision-Making for Future Climate Change (d4PDF) (ref. 26 ). The global climate simulations were conducted using MRI-AGCM3.2 with~60-km grid spacings (ref. 35 ). We conducted largeensemble regional climate simulations that were downscaled from the large ensemble of MRI-AGCM3.2 using NHRCM with 20-km-grid spacing (ref. 36 ). The regional climate simulation covered the East Asian region. The d4PDF's official historical/non-warming datasets covered the period from 1951 to 2010, with calculations extended to the present month as a quasireal-time product.
The HIST simulations were forced by the historical SST and sea ice thickness/concentration based on COBE-SST2 (ref. 37 ) and historical anthropogenic and natural forcing agents such as greenhouse gases and solar irradiance (Representative Concentration Pathway 4.5 emission scenario after 2006) with a 100-member ensemble with different initial conditions and SST perturbations from 1951 to the present. The NonW simulations were forced by historical natural forcing agents and counterfactual "natural" SST and sea ice estimated by removing the warming trends observed in the 20th century. The greenhouse gas concentrations, anthropogenic aerosols and volcanic sulfate aerosols were fixed at 1850 values in the NonW simulation (for more details, refer to ref. 38 ).

Precipitation and circulation indices
The coastal region of Japan's Inland Sea, western Kyushu and eastern Kyushu were defined as coloured grids in Fig. 2a-c, respectively. Depending on the timescale of each event, 72-hour accumulated rainfall (the maximum on each grid during 28 June to 8 July) was used for Case2018, and daily maximum rainfall was used for Case2017 (the maximum during 1-10 July) and Case1993 (the maximum during 1-31 July). Because the rainfall simulated by the 20-km-mesh NHRCM underestimated that which was recorded at the observation site (ref. 18 ), we used probability precipitation with a 50-year return period calculated from the 30-year (1981-2010) d4PDF HIST simulation as a threshold. We estimated the sampling uncertainty in Fig. 4 using the Monte Carlo bootstrap sampling method.
The background circulation index for Case2017 was defined as the difference between the area-averaged 850 hPa geopotential height over the region [10°-30°N, 130°E-160°E] and the region [35°-50°N, 150°E-180°E], which represents a meridional dipole pressure pattern inducing southwesterly moisture flow. The index was calculated every July.