Ultra-High-Resolution Numerical Weather Prediction with a Large Domain Using the K Computer: A Case Study of the Izu Oshima Heavy Rainfall Event on October 15-16, 2013

An intense rainband associated with Typhoon 1326 (Wipha) induced a fatal debris flow on Izu Oshima, Japan, on October 15 – 16, 2013. This rainband formed along a local front between the southeasterly humid warm air around the typhoon and the northeasterly cold air from the Kanto Plain. In this paper, the Japan Meteorological Agency Nonhydrostatic Model was optimized for the “K computer”, and ultra-high-resolution (500 – 250 m grid spacing) numerical simulations of the rainband with a large domain were conducted. Two of main factors that affect a numerical weather prediction (NWP) model, (1) grid spacing and (2) planetary boundary layer (PBL) schemes [Mellor–Yamada–Nakanishi–Niino (MYNN) and Deardorff (DD)], were investigated. Experiments with DD (Exps_DD: grid spacings of 2 km, 500 m, and 250 m) showed better reproducibility of the rainband position than experiments with MYNN (Exps_MYNN: grid spacings of 5 km, 2 km, and 500 m). Exps_DD simulated distinct convective-scale up/downdraft pairs on the southeast/northwest sides of the front, whereas those of Exps_MYNN were not clear. Exps_DD yielded stronger cold pools near the surface than did Exps_MYNN. These differences in the boundary layer structures likely had a large impact on the position of the front and the associated rainband. Exps_DD with the 500-m grid spacing showed the best precipitation performance according to the Fractions Skill Score. To check other factors which influence precipitation forecast, model domain sizes, lateral boundary conditions in nesting simulations, and terrain representations were investigated. In the small domain experiments, the rainCorresponding author: Tsutao Oizumi, 3173-25 Showamachi, Kanazawa-ku, Yokohama, Kanagawa 236-0001, Japan E-mail: oizumi@jamstec.go.jp J-stage Advance Published Date: 30 November 2017 ©The Author(s) 2018. This is an open access article published by the Meteorological Society of Japan under a Creative Commons Attribution 4.0 International (CC BY 4.0) license (http://creativecommons.org/license/by/4.0). Journal of the Meteorological Society of Japan Vol. 96, No. 1 26


Introduction
In Japan, heavy rainfall occasionally causes disasters such as debris flows and floods that can induce severe societal damages. For example, torrential heavy rain caused debris flows on the island of Izu Oshima in October 2013 and in Hiroshima in August 2014. In these two cases, there was a large difference in the spatial scale; the former was a rainband accompanying the stationary front associated with a typhoon, whereas the latter was a line-shaped rainband caused by "back-building formation" in a relatively narrow region. The lengths of the rainbands (exceeding 50 mm 3h −1 at peak period) were approximately 330 and 80 km, respectively. Ogura (1991) investigated local heavy rainfalls in Japan in the 1980s and pointed out that many mesoscale convection systems have a line-shaped precipitation structure. Convective rainbands are normally classified into squall-line and back-building types. In Japan, most torrential rains are caused by backbuilding types in the sense that, viewed from an existing convection cell, new cells are generated in the upstream direction of an environmental wind and the cells develop to form a line-shaped precipitation structure. Tsuguti and Kato (2014) conducted a statistical analysis of synoptic-scale disturbances that cause torrential rains and the shape of precipitation systems in Japan between 1995 and 2009. Their results revealed that the most frequent cases of synoptic-scale disturbances were caused by typhoons and tropical cyclones (32.4 %), followed by stationary fronts, the far-reaching impact of typhoons and tropical cyclones, atmospheric depressions, and cold fronts. Except for the case of typhoons and tropical cyclones, 64.4 % of the precipitation systems were line-shaped. Conversely, Bluestein and Jain (1985) analyzed the organizing processes of mesoscale convective-line systems based on 11-year radar observation data in Oklahoma State and classified the processes into four types: broken line, back building, broken areal, and embedded areal. Seko (2010) also investigated the shapes and main-tenance mechanisms of meso-beta-scale line-shaped precipitation systems in Japan and found their forms to be "side and back-building types", in addition to the squall-line and back-building types.
To mitigate the damages and disasters caused by heavy rain, accurate numerical weather prediction (NWP) with a sufficient lead time is necessary to suggest early evacuation to local residents. The accuracy of NWPs depends on both the initial and boundary conditions and the numerical model. As is well known, the initial condition has a critical impact on the shortrange weather forecast, and data assimilation is widely used to extract information from observations (Kalnay 2003;Lahoz et al. 2010).
Meanwhile, the accuracy of numerical models depends on several factors such as the resolution, domain, dynamics, and physical processes. Finer grid spacing contributes to improving the representation of deep moist convection, reducing finite discretization errors, and expressing more realistic topography. Several studies have focused on the impact of different resolutions in NWPs. For instance, Roberts and Lean (2008) examined whether an increase in the horizontal resolution (grid spacing) of an experiment (12, 4, and 1 km) could produce more skillful precipitation forecasts and concluded that the 1-km experiment had the best performance for rainfall prediction. Weushoff et al. (2010) investigated the performances of three regional NWP models (COSMO, ALADIN, and AROME) over the Alps at different resolutions. Their study showed that all of the high-resolution forecasts outperformed the low-resolution forecasts. Bryan et al. (2003) demonstrated that accurate simulations of convective processes require a horizontal grid spacing on the order of 100 m. Bryan and Morrison (2012) also investigated the sensitivities of microphysics and the horizontal grid distances (4, 1, and 0.25 km) for the simulation of a squall-line and showed that both factors influence the results, even though some quantities (e.g., the surface rainfall) were more sensitive to the resolution.
Many studies have performed sensitivity tests on grid spacing, cumulus parameterization schemes, plan-etary boundary layer (PBL) schemes, and cloud microphysics for various heavy rain events. For example, Tao et al. (2011) examined the sensitivities of several microphysics and PBL schemes in the case of an extreme rainfall event associated with a typhoon in Taiwan using the Weather Research and Forecasting model with 2-km horizontal grid spacing. Their study showed that, the cloud microphysics affected the results rather than the PBL schemes. Fiori et al. (2014) performed a sensitivity experiment on an extreme rainfall event in Italy with a 1-km grid experiment, which used three different convection closures (explicit, Kain-Fritsch, and Betts-Miller-Janjic) and two cloud microphysics schemes. They suggested that the optimal microphysical parameters specified for this event had a positive impact on the prediction of the precipitation intensities. However, they also suggested the necessity to repeat similar experiments for different cases.
In mesoscale modeling, when the turbulence length scale L is comparable to the horizontal resolution Δ, such a region is called the "Terra Incognita" (TI; Wyngaard 2004) or the "gray zone", and it is debatable what type of parameterization is required in such a case. When the horizontal resolution and the boundary layer depth are on similar scales, it is not clear what should be done because both resolved and unresolved turbulent eddies are mixed.
It is well known that closure models for large-eddy simulations (LESs) using schemes such as the Dear dorff (DD; Deardorff 1980) scheme tend to underestimate the subgrid-scale heat flux in the lower boundary layer and overestimate the grid-resolved vertical fluxes of the heat and moisture (the grid-scale disturbance) in cases of coarse-grid spacing (Shin and Hong 2015). In PBL schemes, local schemes that consider local stability only (e.g., the DD scheme) tend to underestimate the vertical mixing and maintain the absolute instability near the surface, which yields overestimations of the grid-scale convection (Ching et al. 2014).
One NWP model, the Japan Meteorological Agency nonhydrostatic model (JMA-NHM; hereafter "NHM"), has implemented the Mellor-Yamada- Nakanishi-Niino (MYNN;Nakanishi and Niino 2004) scheme, which is regarded as one of the most sophisticated PBL schemes available. It is common to use the MYNN scheme when the horizontal resolution is larger than the boundary layer depth (e.g., the horizontal grid spacing is 2 km or more) such as in operational applications. However, when the horizontal resolution becomes sufficiently high to start resolving the boundary layer turbulence, an LES-based model such as the DD scheme seems more suitable. Kato (2011Kato ( , 2012 investigated the dependency of snowfall forecasts on horizontal (0.5 -5 km) and vertical grid spacings and turbulence schemes (MYNN level 3 and DD). In both studies, they found that the turbulent scheme had a larger impact on the precipitation amount than did the horizontal resolution. However, no experiments have been conducted using a high-resolution NWP model to compare MYNN and DD in heavy rain predictions. It is important to examine how MYNN and DD impact the representation of mesoscale convective systems and the quantitative forecasts of intense rain.
Precipitation associated with a frontal system is often enhanced by orographic forcing. Oku et al. (2010) compared for two cloud-resolving models: NHM and the Weather Research and Forecasting model with similar model settings and topography data. They pointed out that not only the differences in model numerics and physics but also slight changes in the representation of topography induce significant differences in the representations of extreme weather. Nunalee et al. (2015) used a high-resolution numerical model (1-km grid spacing) to simulate actual cases of atmospheric flow past mountainous islands using three different topography datasets. They found that the use of accurate orographic boundary conditions is more advantageous when running high-resolution mesoscale models. Roberts et al. (2009) investigated the orographic enhancement of rain in England using the Met Office Unified Model (UM) with different grid spacings of 12, 4, and 1 km. Their results showed the benefit of increased resolution in the UM and the benefit of coupling high-resolution rainfall forecasts to the probability distributed model for flood warnings.
In general, the model domain size is determined by a trade-off between the spatiotemporal scales of the targeted phenomena and the available computational resources. Several studies have investigated the influence of the domain size on precipitation predictions in climate models. Larsen et al. (2013) examined the regional climate model (HIRHAM) with different model domain sizes (1400 -5500 km) and different horizontal grid spacings (5.5, 11, and 12 km) for precipitation simulations. They concluded that, the domain sizes were more important than the resolutions in reproducing the precipitation. Bray et al. (2011) investigated the sensitivity of the domain size and the buffer zone to the uncertainty of extreme rainfall downscaling using the Pennsylvania State University-National Center for Atmospheric Research fifth-generation mesoscale model (MM5). They noted that the domain size and the buffer zone have a significant impact on the model rainfall estimates. Davies (2014) determined some general guidelines for lateral boundary condition (LBC) strategies. These strategies imply that, for a small domain, forecast fields are largely determined by LBCs and as a consequence low update frequencies of LBCs cause a loss of information in forecast fields. Wang et al. (2016) tested the sensitivity of heavy precipitation to various model configurations using the Australian Community Climate and Earth-System Simulator via a series of convection-permitting simulations (0.5 -4 km) for three heavy precipitation cases. They showed that impacts of model resolutions and domain sizes on quantitative precipitation forecast can be compared with impacts of physical parameterization schemes and initial conditions. These previous studies independently suggest the benefits of finer resolutions and model domains in selected cases. However, no study has been conducted for ultra-high-resolution simulations (100-m scale) with large model domains for NWP. Such a high-resolution simulation with a large domain experiment requires huge computational resources. However, using a large domain has the benefit of reducing the impact of the boundary conditions.
In Japan, the "K computer" owned by RIKEN, whose performance in the LINPACK (LINear equations software PACKage) benchmark floating-point operation was ranked first in the world's top supercomputers in June 2011, and has been used for research since September 2012. Half of the computational resources of the K computer has been assigned to "High-Performance Computing Infrastructure Strategic Program for Innovative Research (SPIRE)" funded by the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan. As a sub subject of one of the five strategic fields of SPIRE, "Ultra-high Precision Meso-Scale Weather Prediction" has been conducted  with three subthemes: 1) the development of cloud-resolving fourdimensional data assimilation systems, 2) the development and validation of a cloud-resolving ensemble analysis and forecast system, and 3) basic research using very high-resolution atmospheric models. Several of the published results (e.g., Kobayashi et al. 2013;Kunii 2014a, b;Chen et al. 2015a, b;Duc et al. 2015;Ito, K. et al. 2015;Seko et al. 2015;Ito and Niino 2016) are on the project website (https://www. jamstec.go.jp/hpci-sp/en/strategy/mswp.html).
The objectives of this study are to conduct an ultra-high-resolution experiment (500 to 250-m grid spacing) with a large domain and to examine important factors for accurate forecasts of heavy rain events. The target event is the Izu Oshima heavy rain that occurred in October 2013 when a strong typhoon approached East Japan. The heavy rain was brought on by a rainband associated with the stationary front, which was maintained by the generation of new convection cells on the inflow (warm) side; however, its horizontal scale of approximately 330 km was larger than those described in previous studies of line-shaped precipitation with back-building formation.
To accomplish ultra-high-resolution experiments with a large domain, first we optimized NHM for the K computer and modified the pre-and post-processes of the NHM.
The important factors, (1) the grid spacing (5 km, 2 km, 500 m, and 250 m) and (2) the PBL scheme (MYNN levels 2.5 and 3 and DD), for heavy rain prediction were investigated using the optimized NHM with a large domain (1600 km × 1100 km). To discuss the large model domain impacts, we investigated the impacts of experiments with a small model domain (200 km 2 ) and LBCs on the rainband formation. The sensitivity analyses of the model terrain were conducted with two terrain datasets made from different digital elevation maps. To obtain rigorous results, some significant factors were examined with respect to different initial times . The obtained results were quantitatively evaluated using the Fractions Skill Score.
This paper is organized as follows. In Section 2, the details of the Izu Oshima heavy rain event in October 2013 are described. In Section 3, the optimization of the NHM for the K computer is described, and the design of the experiments is presented. Section 4 shows the results of the simulations, and the main factors determining the rainband position are examined. In Section 5, the impacts of the model domain size and the LBCs and the different representations of the terrain are discussed. The robustness of this study was investigated for different initial times, and a qualitative verification of the precipitation is shown. The summary and concluding remarks are given in Section 6.

The heavy rainfall on Izu Oshima Island on October 16, 2013
Izu Oshima is a volcanic island located approximately 120 km south of the central Tokyo Metropolis area and approximately 30 km east of the Izu Peninsula. In the early dawn of JST (Japan Standard Time, + 9 GMT) October 16, 2013, Typhoon 1326 (Wipha) approached Izu Oshima, and torrential heavy rain caused massive landslides on the island. The damage statistics reported that there were 39 dead and/or missing persons and 177 affected households.
The typhoon formed near the Mariana Islands at 2100 JST on October 10, 2013. The typhoon initially moved westward and then curved to the northeast on the morning of October 15. In the early morning hours of October 16, the typhoon was closest to Izu Oshima, and the rainband associated with the typhoon brought very heavy rain over Izu Oshima. After that, the typhoon moved northeastward. Figure 1 shows the synoptic weather chart at 0300 JST on October 16. The center of the typhoon was located approximately at a latitude of 32.8°N and a longitude of 138.6°E, and the minimum sea-level pressure was 950 hPa. The synoptic chart depicts a stationary front extending from southwest to northeast at 2100 JST on October 15 (JMA 2014); then, the front slowly shifted northward, corresponding to the northeastward movement of the typhoon.
The Japan Meteorological Agency (JMA 2014) analyzed the mechanism of the heavy rain event on Izu Oshima. A stationary front formed between the humid warm air from the south supplied by the typhoon and the northeasterly cold air from the Kanto Plain. An intense rainband formed and was maintained along the stationary front. Kato (2015) classified the rainband as a "broken-line type"; however, Bluestein and Jain (1985) and Ogura (1997) stated that broken-line-type formations tend to occur along cold fronts. In the present case, the precipitation band occurred along the stationary front associated with the typhoon. Therefore, the situation surrounding the occurrence of the Izu Oshima heavy rain was different from those in the previously mentioned studies. Hara (2014) analyzed the MSM forecast results for the Izu Oshima heavy rain event and noted that a coastal front formed between the cold air over the Kanto Plain and the humid warm air from the typhoon. Bosart and Carr (1978) and Bosart and Dean (1991) investigated the formation mechanism of a coastal front in the eastern United States ahead of a tropical cyclone and noted the importance of surface convergence driven by land-sea thermal contrasts. Figure 2a shows the observed one-hour precipitation from 2200 JST October 15 to 0600 JST October 16 (JMA's precipitation analysis by Radar-AMeDAS (Automated Meteorological Data Acquisition System of JMA)) and other rain gauge observations (JMA 2013). Note that the accuracy of Radar-AMeDAS data on the ocean area is lower than that of the data on the land area because it is far from meteorological radars and no rain gauge observations are available over the ocean. The rainfall around Izu Oshima started at 0900 JST on October 15 and gradually increased in intensity. As shown in Fig. 2a, a rainband formed over Izu Oshima around 2300 JST on October 15 and then intensified. The highest period of rain over Izu Oshima was from 0100 JST to 0500 JST on October 16. After 0500 JST, the rainband weakened. Figure 2b shows the three-hour accumulated precipitation (0100 -0400 JST on October 16) during the peak period of the event. Izu Oshima was located at the center of the intense rainband, which extended from the south of the Izu Peninsula to the Boso Peninsula.
As seen in the enlarged view (Fig. 2c), the threehour accumulated precipitation on Izu Oshima was more than 200 mm on the northern half of the island and 150 -200 mm on the southern half. The most intense rain was analyzed on the northeast part of the island (the point marked by the white cross). Figure 3 shows the topography of Izu Oshima. The shape of the island is nearly oval and elongated in the NNW-SSE direction, with a major axis of approximately 15 km and a minor axis of approximately 9 km. Mt. Mihara (whose highest peak is 764 m above sea level) is located at the central part of the island.
We targeted the heavy rain event on Izu Oshima from 1800 JST October 15 to 0600 JST October 16, 2013. An AMeDAS station at Oshima (P1 in Fig.  3) recorded a three-hour accumulated precipitation of 328 mm during the peak period of the event (0100 -0400 JST) and a 24-hour precipitation of 824 mm at 0800 JST on October 16. The recorded 24-hour accumulated precipitation was 2.5 times larger than the average monthly precipitation at the same observation point. Another AMeDAS station at Kitanoyama (P2 in Fig. 3) recorded a three-hour accumulated precipitation of 170 mm at the peak period 1 .

Optimization of the NHM for the K computer
In this study, the NHM (Saito et al. 2006(Saito et al. , 2007Saito 2012) for the K computer was optimized to conduct these enormous computations. The K computer is a scalar-type computing system with 88,128 nodes (82,944 computational nodes and 5184 I/O nodes). NHM was originally developed for a vector-type parallel computing system at JMA; therefore, tuning its performance was critical to efficiently run the highperformance computing system.
The authors and Fujitsu Co. Ltd modified the NHM's FOTRAN90 source codes for the K computer. First, we monitored the computational cost of each computational loop and detected 160 expensive loops, whose costs comprised approximately 77 % of the total computational time. The status of the thread parallelization and the single instruction multiple data process was analyzed for each loop. Then, we selected the loops that could be potentially improved, and five optimization techniques were applied to the 51 loops. A more detailed description of this optimization process is given in Oizumi et al. (2015). Figure 4 compares the peak performance and the elapsed time between the original NHM and the optimized NHM (hereafter, "KNHM"). The measurement was carried out using 72 and 288 nodes of the K computer. The measurement conditions were as follows: CASE 1 measured only the time integration part in the main loop, whereas CASE 2 measured the entire elapsed time of the main loop in the actual computation of the two experiments. The K computer compiler, which was released in February 2014, was employed.
Using KNHM, the peak performance of CASE 1  with 72 nodes (Fig. 4a) was improved from 4.6 % in the original NHM to 5.5 %, and the elapsed time ( Fig. 4e) was reduced from 181 to 159 s. The peak performance of CASE 1 with 288 nodes (Fig. 4b) was improved from 4.1 % to 4.6 %, and the elapsed time ( Fig. 4f ) was reduced from 52 to 48 s. A similar improvement was obtained for CASE 2 with 72 and 288 nodes in both the peak performance (Figs. 4c, d) and the elapsed time (Figs. 4g,h). CASE 2 shows less improvement than CASE 1 because the elapsed time of the input-output processes was included. The peak performance improvement with 72 nodes was more distinct than that with 288 nodes because the latency for the communications between 288 nodes is larger than that between 72 nodes for the same model domain. The weak scalability (a scalability test fixes the model size per computational node) of KNHM was 97 %. We tested KNHM with all 82,944 nodes of the K computer and confirmed that the output data from 72 and 82,944 nodes are identical at the binary level.
The pre-and post-processes and file IO of the NHM were also modified to effectively use the K computer (Oizumi et al. 2015).

Design of the experiments
KNHM was used as the NWP model to simulate the heavy rain event on Izu Oshima. Figure 5 shows the two different model domains (referred to as LARGE and SMALL). Here LARGE covers a domain of 1600 km × 1100 km from the south part of the Tohoku district to Kyushu, which is the same as that of the JMA's Local Forecast Model (LFM: JMA 2013) previously operated from October 2010 to May 2011. Experiments with the SMALL domain (a 200-km 2 region whose center is Izu Oshima) were conducted to determine the impact of the model domain size and the LBCs on the simulation results. The line AB is the position of the vertical cross sections of the cloud water, rainwater, and potential temperature shown later. The 20-km2 regions C and D indicate the sampling areas of the vertical profiles in Section 4.2. The square E is the precipitation verification domain for the Fractions Skill Score.

a. Experimental settings for the LARGE domain
We tested three initial times at 1800 JST and 2100 JST on October 15 and at 0000 JST on October 16; however, most experiments were started at 2100 JST on October 15. The forecast lengths were 12, 9, and 6 hours, respectively (until 0600 JST on October 16). The mesoscale analyses of JMA (MA; JMA 2013) were used as the initial conditions. These analyses were also used for the three-hourly LBCs to reduce the uncertainty in the simulation due to the quality of the lateral boundary data.
The details of the model configuration and the experiment names are given in Table 1. The 5-km grid experiment used nearly the same settings as the JMA operational Mesoscale Model (MSM: 5-km grid spacing, 50 vertical layers, and a model top of 21,800 m). The 2-km grid experiment used the same settings as the LFM of April 2014 (LFM: 2-km grid spacing, 60 vertical layers, and a model top of 20,189 m). The 500-and 250-m grid experiments used the same model top as the LFM, and the number of vertical layers was increased to 85 and 168, respectively.
The physics processes of the 5-km grid experiment followed the MSM settings, and those of the other grid experiments followed the LFM settings. The model included the bulk cloud microphysics introduced by Ikawa et al. (1991). The scheme predicts the mixing ratios of six water species (water vapor, cloud water, rain, cloud ice, snow, and graupel) and the number concentration of the cloud ice. In addition, the 5-km grid experiment used a modified version of the Kain-Fritsch cumulus parameterization scheme (Kain and Fritsch 1993).
Three turbulence closure schemes, Mellor-Yamada-Nakanishi-Niino level-3 (MYNN) and level-2.5 (MYNN25) Niino 2004, 2006) and Deardorff (DD: Deardorff 1980), were employed in this study. JMA uses MYNN for LFM for operation. In MSM, MYNN was used for operation until the early 2015. After that, MYNN25 has been employed to improve the forecast of winter precipitation. In this study, most experiments were conducted with MYNN, and in addition, we tested MYNN25 for the 5-km grid experiment. To investigate the impact of the model grid spacing and the turbulence closure models (MYNN and DD) on the simulation results, the 2-km and 500-m grid experiments were examined. The 250-m grid experiments used only DD, assuming that the horizontal grid spacing was smaller to the boundary layer depth.

b. Experimental settings with the SMALL domain
The impact of the model domain size and the LBCs on the simulation results was investigated using the 2-km and 500-m grid experiments. Table 2 shows the experiment names and the model settings. All model settings of the 2-km and 500-m grid experiments were the same as those of CM2kmMY and CM500mDD with LARGE, respectively, except that the model domain was SMALL. The experiment names with "_small" were the same as the experiments with the large domain except for the model domain size, and these performed nine-hour simulations starting at 2100 JST on October 15. All experiment names with "_nest" used a one-way nesting procedure and performed six-hour simulations starting at 0000 JST on October 16. CM2kmMY_nest and CM500mDD_ nest were one-way nested inside CM5kmMY and CM2kmMY, respectively. In the nesting procedure, the LBCs of the nested models were given by the spatiotemporal interpolation of the parent model forecasts with one-hour updates. The experiments CM2kmMY_ nest_noBND and CM500mDD_nest_noBND were the same as CM2kmMY_nest and CM500mDD_nest, respectively, except that the cloud microphysical quantities were not included in the LBCs. The experiments CM2kmMY_nest_3hr and CM500mDD_nest_3hr were the same as CM2kmMY_nest and CM500mDD_ nest, respectively, except that the update frequency of LBCs was three hours.

c. Terrain data
The terrain data of MSM was generated via an interpolation of the GTOPO30 (Global 30 arc-second elevation with a horizontal grid spacing of approximately 1 km) dataset. In the operational forecast, to avoid the risk of numerical problems, the maximum model terrain inclination between two adjacent grid points was limited to 15 % by a control parameter (GRMAX) in the preprocessing when preparing the model terrain. Figure 6 shows the model terrain of Izu Oshima with various grid resolutions: Figs. 6a, 6b, 6c, and 6d illustrate the terrains of the 5-km, 2-km, 500-m, and 250-m grid experiments, respectively, interpolated from GTOPO30 with a GRMAX of 15 %, and Figs. 6e and 6f show the terrain of the 500-and 250-m grid experiments, respectively, interpolated from KTOPO (a digital elevation model with 50-m grid spacing provided by the Geospatial Information Authority of Japan) without a limit on the inclination. In Fig.  6a, the 5-km grid experiment terrain is considerably smoother than the actual terrain (Fig. 3). In Fig. 6b, the highest peak of the 2-km grid experiment (497 m) is approximately 260 m lower than that of the actual terrain. In the 500- (Fig. 6c) and 250-m (Fig. 6d) grid experiments, there is no significant difference between the terrains because the grid spacing of GTOPO30 is approximately 1 km. The contour shapes of the 500- (Fig. 6e) and 250-m (Fig. 6f ) grid experiments represent the characteristics of the actual terrain well, as shown in Fig. 3; in particular, Fig. 6f shows the most realistic terrain with the highest peak of 693 m.
The experiment name CM500mDD_no_Izu indicates an experiment where the terrain of Izu Oshima is removed from the LARGE domain. CM500mDD_K and CM500mDD_K used terrain data that interpolated KTOPO without a limit on the inclination. The other model settings are the same as those in Table 1. Figure 7 shows the mean sea-level pressure at 0300 JST on October 16 according to the experiments with the large model domain started at 2100 JST on October 15 (Table 1). The predicted center positions of the typhoon are located within 32.5 -33.0°N and 138.0 -138.5°E, which agree well with the JMA best track data (32.8°N, 138.6°E). The range of the central pressures of the typhoon was 956 -959 hPa, which was 6 -9 hPa higher than the central pressure of the best track data (950 hPa). The experiments with DD tended to predict slightly lower central pressures than the experiments with MYNN. Figure 8 shows the three-hour precipitation (0100 -0400 JST, October 16) of the JMA precipitation analysis and the corresponding results of the numerical models. In the precipitation analysis (Fig. 8a, enlarged view of Fig. 2a), Izu Oshima is embedded at the center of an intense rainband, which extends from the Boso Peninsula to the south of the Izu Peninsula. As previously shown in Fig. 2b, the precipitation intensity was 150 -200 mm on the southern part of the island and more than 200 mm on the northern part of the island.

The location of the typhoon and precipitation distributions
In all experiments with the MYNN scheme (KF 5kmMY) (Fig. 8b), KF5kmMY25 (figure not shown, nearly the same as that in Fig. 8b), CM2kmMY (Fig. 8c), and CM500mMY (Fig. 8e), the rainband position was predicted further northwest compared to the observations (Fig. 8a). In KF5kmMY (Fig. 8b) and KF5kmMY25, the precipitation intensity over Izu Oshima was 70 -100 mm on the northern part of the island and 50 -70 mm on the southern part. In CM2kmMY (Fig. 8c) and CM500mMY (Fig. 8e), the precipitation over the island was 70 -120 mm on the northern part and 50 -70 mm on the southern part.
Conversely, in the experiments with DD (CM2km DD) (Fig. 8d), CM500mDD (Fig. 8f ), and CM250m DD (Fig. 8g), the predicted rainband covered Izu Oshima as in the observations. CM2kmDD (Fig. 8d) reproduced the intense rain (150 -200 mm) over the center of the island and the 100 -150 mm precipitation over other areas. In CM500mDD (Fig. 8f ), the precipitation intensity was 120 -240 mm on the northern half of the island and 70 -100 mm on the southern half. CM250mDD (Fig. 8g) showed a similar precip- itation distribution on the island to Fig. 8f, where the precipitation was 150 -230 mm on the northern part of the island and 70 -100 mm on the southern part. The rainband in the experiments with the MYNN scheme deviated northwest from the observations regardless of the grid spacing, and the position of the rainband was considerably improved in the experiments with DD. This result indicates that the turbulence closure model has a considerable impact on rainband formation. We investigate the causes of this difference and it is explained in the next subsection.
CM500mDD and CM250mDD reproduced the peak of the intense rain better than CM2kmDD; however, the areas of intense rain and the rainband position were seemingly not very different. In Section 5.4, we show a quantitative evaluation of the simulated precipitation using the Fractions Skill Score.

The rainband formation
The meteorological situation of the Izu Oshima heavy rainfall event was promptly reported by JMA (2014); the intense rainband formed along the front between the southeasterly warm and moist airflow around Typhoon 1326 (Whipa) and the northerly cold airflow from the southern part of the Kanto Plain. In this section, the factors significant to the formation and maintenance of the front and the rainband are investigated with the experiments CM2kmMY, CM2kmDD, CM500mMY, and CM500kmDD. We compared the near surface temperature and the wind distribution, the vertical flow distribution at the low level (altitude 1 km), the mixing ratios of the cloud water and rainwater, and the potential temperature.
a. Horizontal and vertical distribution of the model quantities Figure 9 shows the horizontal distribution of the temperature and the horizontal wind at an altitude of 20 m (the lowest level of the model). The upper four panels show the results of CM2kmMY (Fig. 9a1), CM2kmDD (Fig. 9b1), CM500mMY (Fig. 9c1), and CM500mDD (Fig. 9d1) at 0100 JST on October 16, just prior to the intense rain on Izu Oshima, and the lower four panels (Figs. 9a2 -d2) show the results at  0300 JST on October 16, when the rainfall was the most intense. In all figures, the overall trends of the temperature distributions and wind directions were similar. The temperature of the southeasterly flow south of the front was 24 -27°C, whereas the northeasterly air from the Kanto Plain north of the front was below 21°C. The front extending from northeast to southwest formed between these two flows with large differences in temperature and distinct horizontal shears. The wind directions of the warm and moist airflows were nearly perpendicular to the front, whereas on the cold air side, the wind directions were northeasterly, and nearly parallel to the front.
In CM2kmMY, the front overlapped with Izu Oshima at 0100 JST (Fig. 9a1); then, the front shifted westward and was located between Izu Oshima and the Izu Peninsula at 0300 JST (Fig. 9a2). CM500mMY (Figs. 9c1,c2) show similar results to CM2kmMY for the position of the front and the temperature distribution. Meanwhile, in CM2kmDD, the front was simulated southeast of the island at 0100 JST (Fig. 9b1), and the front subsequently remained southeast of Izu Oshima (Fig. 9b2). In CM500mDD (Figs. 9d1, d2), the positions of the front were slightly closer to the island than those of CM2kmDD. Note that, in the experiments with DD, the temperatures west of the Tokyo Bay and over the western part of the Boso Peninsula were 1 -2°C lower than those in experiments with MYNN. These results indicate that the PBL scheme significantly affected the position of the front.
The upper four panels of Fig. 10 show the horizontal distributions of the vertical flow at an altitude of 1 km and the vertical cross section of the cloud water and rainwater along the line AB (the same line AB as in Fig. 5) at 0100 JST on October 16 as modeled by CM2kmMY (Figs. 10a1, a2) and CM2kmDD (Figs.  10b1, b2). The arrow symbols indicate the position of the simulated front. In CM2kmDD (Fig. 10b1), a distinct band-shaped updraft is seen along the front, whereas in CM2kmMY (Fig. 10a1), no apparent band-shaped updraft was simulated along the front. In CM2kmDD, the updraft along the front appeared in an early stage of the simulation, i.e., at 2130 JST (figure not shown).
In the vertical cross section, the cloud water in CM2kmDD (Fig. 10b2) developed along the updraft on the southeast side of the front (arrow), and the rainwater appeared on the northwest side of the front. In CM2kmMY (Fig. 10a2), the cloud water appeared from the more southeastern side but did not show any concentration at a particular position, and only weak rains appeared along the front.
With the finer grid spacing (Figs. 10c1 -d2) in both CM500mMY and CM500mDD, the vertical flows were stronger and appeared in smaller cells. In CM 500mDD (Fig. 10d1), convective-scale updrafts and downdrafts appeared southeast and northwest of the front, respectively. In CM500mMY (Fig. 10c1), however, the updraft along the front was weak, and the downdraft was obscure.
In CM500mDD (Fig. 10d2), the cloud water appeared in a more southeastern position than in CM 2kmDD; however, a large amount of the cloud water was simulated at the position of the strong updraft along the front, and the rainwater was associated with the downdraft on the slightly northwest side of the front. These characteristics were not clear in CM-500mMY (Fig. 10c2). Figure 11 is the same as Fig. 10 except for the valid time which is 0300 JST on October 16than that of CM500mDD (the dashed red line) below/above 500 m. The vertical gradients of the potential temperature in the range between 300 m and 2.0 km of CM500mMY were smaller than those of CM500mDD, which indicates that the profile of CM500mMY was more upright at this height. This tendency of the vertical gradients is similar to that of square C; however, the difference between the two schemes was larger in square D than in square C. One possible reason for this is the differences in the temperatures between the sea surface and the lowest atmosphere in the two squares. The sea surface potential temperature is 298.8 K for square C (indicated by an arrow symbol in Fig.  13a) and 298.4 K for square D. Therefore, the difference in the temperatures between the sea surface and the lowest atmosphere in square C is small, whereas in square D, the lowest atmosphere is absolutely unstable owing to the heating from the relatively warm sea surface. Sensible heat flux at the sea surface is even less than zero in square C, and in square D, the flux with DD had larger values than that with MYNN (see Supplement 1). These results indicate that vertical mixing in the MYNN experiments tends to be larger than that in DD for an unstable atmosphere, which likely yields the larger difference in the vertical profiles of the (potential) temperature compared to that found in square C. Another possible reason is the differences in the positions of the simulated front, which was located in a more northwest position in CM500mMY than in CM500mDD; therefore, the air in square D was likely affected more by the northeasterly flow from the Kanto Plain in CM500mDD. The weak  cold pool formed in the northwest side of the front in CM500mMY (see Fig. 12c) might also contribute to the larger differences in the vertical profiles of the potential temperatures of the two schemes in square D compared to those in square C. Figure 13b shows the vertical profile of the mixing ratio of the water vapor (qv). In square C, the qv of CM500mDD (the bold red line) was slightly larger than that of CM500mMY (the bold black line) below 500 m, and the qv of CM500mMY was larger than that of CM500mDD above 500 m. The vertical gradient of qv in CM500mMY was slightly smaller than that in CM500mDD at the altitude between 300 m and 2.0 km, which is a similar tendency to that of the potential temperature (Fig. 13a). In square D, the vertical profile of the qv of CM500mMY (the dashed black line) was more upright than that of CM500mDD (the dashed red line).
In Fig. 13c, the vertical profiles of the equivalent potential temperatures (θe) are shown. In square C, the difference in θe values between CM500mMY and CM500mDD is small because the tendencies of the potential temperature and qv are canceled out The θe of CM500mDD (the bold red line) is slightly higher than that of CM500mMY (the bold black line) below 500 m, and above 500 m, the θe of CM500mMY is slightly higher than that of CM500mDD. The vertical profile of CM500mMY is slightly more upright than that of CM500mDD, which indicates that the atmosphere in CM500mDD is more unstable. In square D, θe was positively affected by the tendency of the magnitude of the potential temperature and qv. The θe of CM500mMY (the dashed black line) was larger than that of CM500mDD (the dashed red line) below 500 m, and the θe of CM500mDD was larger than that of CM500mMY in the altitudes between 500 m and 1.7 km. The above results indicate that the MY scheme tends to evaluate the vertical mixing in the boundary layer as being larger than DD does and that the thermal instability in the lower boundary layer is maintained more in DD, which is consistent with previous studies (Ching et al. 2014;Shin and Hong 2015). Figure 13d shows the vertical profiles of the absolute value of the horizontal wind speeds [(u 2 + v 2 ) 1/2 ]. As described in Fig. 9, the southeasterly wind direction in square C (the warm and moist side) is nearly perpendicular to the front, whereas in square D (the cold air side), the northeasterly wind direction is nearly parallel to the front. In square C, the wind speed of CM500mMY (the bold black line) was smaller than that of CM500mDD (the bold red line) between 200 m and 1.5 km, and the vertical gradient of CM500mMY was smaller in the lower part of the mixing layer below 600 m. In square D, the wind speed of CM500mDD (the bold red line) was generally larger than that of CM500mMY (the bold black line). Figure 14 shows the distribution of the equivalent potential temperature at the lowest model level (20 m above ground level) and at an approximate altitude of 2 km at 0100 JST on October 16. In square C (the warm and moist side), the difference in θe values between CM500mMY and CM500mDD was small at both levels. Conversely, in square D (the cold air side), the θe of CM500mDD was lower than that of CM500mMY. Square D in CM500mDD was more distant from the front than in CM500mMY; therefore, the air in square D is likely more affected by the cold air flow from the Kanto Plain.
In Supplement 2, we show the vertical profiles of θ, qv, θe, and the wind speed in the cold air at a different position. In this figure, the profiles CM500mDD_D' are computed in the square D' (e) so that the distance from the front has the same magnitude as that for CM-500mMY_D in Fig. 13. The tendencies of the vertical profiles of CM500mDD_D' are similar to those of CM500mDD_D in Fig. 13.
In this case study, the experiment with DD showed better performance than the experiment with MYNN. However, these results do not indicate that the DD scheme is more appropriate than the MYNN scheme. In general, the direct application of the Deardorff scheme is only justified for models with LES resolutions. In particular, it is necessary to use the 2-km grid spacing model with DD.

Impact of the model domain size and the LBC
In Section 4, the finer grid experiments were better able to reproduce the observed rainband in terms of its position and rainfall intensity over Izu Oshima. However, a very high-resolution experiment requires huge computational resources such as the K computer. For this reason, most previous studies have conducted high-resolution experiments with nesting procedures in limited model domains. In this subsection, we examine the impact of the model domain size and the LBC on the simulation of heavy rain. For comparison, the SMALL domain shown in Fig. 5 was used with the model settings of CM2kmMY and CM500mDD. In the experiments with the small model domain, the effect of the one-way nesting procedure was also examined with and without the inclusion of the cloud microphysical quantities in the LBC. The names of the experiments and their specifications are listed in Table 2. Figure 15 shows the three-hour accumulated precipitation for 0100 -0400 JST on October 16 for the experiments with a small domain. The upper four panels show the results of the 2-km experiments, and the lower four panels show those of the 500-m experiments corresponding to CM2kmMY (Fig. 8c) and CM500mDD (Fig. 8f ), respectively.
In CM2kmMY_small (Fig. 15a) and CM500mDD_ small (Fig. 15e) without nesting, the precipitation patterns are very different from those of the experiments with a large domain (CM2kmMY and CM500mDD). The intense rainband does not appear between the Izu Peninsula and the Boso Peninsula, and a different unrealistic isolated intense rain area appears around the Tokyo Bay and the Boso Peninsula.
In CM2kmMY_nest (Fig. 15b), the features of the rain distribution around Izu Oshima were similar to those in CM2kmMY even though the overall shape of the rainband was slightly different. The rainband in CM500mDD_nest (Fig. 15f ) also showed similar features to the one in CM500mDD.
In the experiments where the cloud microphysical quantities were excluded from the LBCs (CM2km-MY_nest_noBND (Fig. 15c) and CM500mDD_nest_ noBND (Fig. 15g), the overall shapes of the rainband were similar to those of CM2kmMY_nest (Fig. 15b) and CM500mDD_nest (Fig. 15f ), respectively. However, the precipitation near the inflow boundary (the lower right of the figures) of the model domain was significantly decreased, whereas the precipitation in the southwestern part of the rainbands was increased in both experiments.
In CM2kmMY_nest_3hr (Fig. 15d) and CM500m DD_nest_3hr (Fig. 15h), the impacts of the update frequency of the boundary conditions were checked by changing the time interval of the update from one hour to three hours. The intensity of the rainband of CM2kmMY_nest_3hr was weakened compared to that of CM2kmMY_nest (Fig. 15b). In CM500mDD_ nest_3hr (Fig. 15h), the rainband became indistinct, and a fake intense precipitation similar to that in CM 500mDD_small appeared around the Boso Peninsula.
For comparison, additional experiments were conducted with the CM2kmMY setting using a different size domain, a 400-km 2 region centered on Izu Oshima (Supplement 3). In all experiments, the shape of the rainband was improved compared to that of the experiments with the SMALL domain, particularly in CM2kmMY3_small_ext and CM2kmMY_nest_3hr_ ext. A reduction in the precipitation near the inflow boundary is still seen in CM2kmMY_nest_noBND_ ext; however, it does not seriously affect the features of the rainband. The small domain model is more sensitive to the LBC, and the inclusion of the cloud microphysical quantities in the LBC is critical if the model domain is smaller than the target phenomenon scale. The good similarity between the original simulations (Figs. 8c, e) and the nested simulations (Figs. 15b, f ) suggests that the nesting procedure of the NHM works properly and that the basic features of the event are well reproduced with the small model domain if the cloud microphysical quantities are included and the boundary is updated hourly. However, the required model domain size should depend on the spatiotemporal scale of the event.

Impact of terrain on the precipitation distribution
One plausible cause of the intense rain on Izu Oshima is the orographic effect of the island. In this subsection, the impact of the terrain on the precipitation distribution is investigated for the experiments CM500mDD and CM250mDD, which showed relatively good performance in this study. First, the impact of the terrain of Izu Oshima on the rainband position and precipitation intensity is examined for CM500mDD with the LARGE domain.

a. Effect of terrain of Izu Oshima on the rainband
The experiment CM500mDD_no_Izu used orographic data where the terrain height of Izu Oshima was set to zero with the LARGE domain. Here, the distribution of land and sea was maintained as in the original model setting (Fig. 6c). Figure 16 shows the three-hour accumulated precipitation from 0100 JST to 0400 JST on October 16 by CM500mDD_no_Izu. The overall features of the rainband (the position and precipitation area) in this experiment were nearly the same as those in CM500m DD (Fig. 8f ); however, the precipitation intensity around the island was very different. The precipitation intensity in the north and central parts of the island in CM500mDD_no_Izu was reduced from over 200 mm to less than 120 mm. These results clearly indicate that the precipitation on Izu Oshima was locally enhanced by the orographic effects of the island.

b. Experiments with finer terrain data
In this study, the horizontal grid spacing of the highest resolution experiment was 250 m. However, as we mentioned in Section 3.2c, the terrain data were given by interpolation from GTOPO30 (0.1° grid spacing), and the upper limit of the terrain slope was set to 15 % (GRMAX = 0.15); therefore, the 250-m grid spacing model terrain (Fig. 6d) was nearly the same as that used in CM500mDD (Fig. 6c). One possible reason that the three-hour accumulated precipitations of CM500mDD (Fig. 8f ) and CM250mDD (Fig.  8g) showed a small difference was likely to be the terrain representation. In this subsection, we examine the impact of the finer terrain data interpolated from KTOPO (50-m grid spacing digital elevation model without a limit on the inclination) for the experiments with 500-and 250-m grid distances. Figure 17 shows the three-hour accumulated precipitation from 0100 JST to 0400 JST on October 16 around Izu Oshima. Here, CM500mDD (Fig. 17a) and CM250mDD (Fig. 17c) are enlarged views of Figs. 8f and 8g, respectively, and CM500mDD_K (Fig. 17b) and CM250mDD_K (Fig. 17d) used the terrain data from KTOPO (Figs. 6e, f, respectively). The intense rain areas (orange: 120 -150 mm, red: 150 -200 mm, and pink: more than 200 mm) of CM500mDD_K and CM250mDD_K were spread out more than those of CM500mDD and CM250mDD, respectively. In particular, CM250mDD_K simulated precipitation intensity similar to the analyzed rain (Fig. 2b), where precipitation of more than 200 mm covered the northern half of the island and the most intense rain (277 mm) was on the northeast part of the island. These results indicate that the precipitation intensity over Izu Oshima is better reproduced when a finer grid spacing and terrain representation of the island is used.

Experiments with different initial times
The experiments in the previous sections showed the results of nine-hour forecasts starting at 2100 JST on October 15. To confirm the robustness of the  tendency obtained by the experiments in the previous sections, a six-hour forecast starting at 0000 JST on October 16 and a 12-hour forecast starting at 1800 JST on October 15 were additionally carried out and the results are analyzed in this subsection. The experimental setting is the same as that in CM2kmMY, CM2kmDD, CM500mMY, and CM500mDD listed in Table 1. Figure 18 shows the three-hour accumulated precipitation from 0100 JST to 0400 JST on October 16 corresponding to Fig. 8. The upper panels show the results of the six-hour forecasts starting at 0000 JST on October 16, whereas the lower panels show the 12hour forecasts starting at 1800 JST on October 15. In the six-hour forecasts, the positions of the rainband in CM2kmDD (Fig. 18b) and CM500mDD (Fig. 18d) were very similar to those in CM2kmDD (Fig. 8d) and CM500mDD (Fig. 8f ), respectively, and more intense rain appeared over the island. In CM2kmMY_6hr (Fig.  18a) and CM500mMY_6hr (Fig. 18c), the rainband positions were slightly closer to Izu Oshima than those of CM2kmMY (Fig. 8c) and CM500mMY (Fig.  8e), respectively, but still northwest of Izu Oshima.
In the 12-hour simulations using MYNN, the rainband was simulated further northwest, as shown by CM2kmMY_12hr (Fig. 18e) and CM500kmMY_12hr (Fig. 18g). In these experiments, intense rain appeared only over the Izu Peninsula as an isolated area. Conversely, the position of the rainband in CM2km DD_12hr (Fig. 18f ) was nearly the same as that in CM2kmDD (Fig. 8d). In CM500mDD_12hr (Fig.  18h), the rainband position was slightly shifted to the northwest, but intense rain over the island was still simulated.
In the above sensitivity test, experiments with the MYNN scheme showed a larger dependency on the initial times than the experiments with the DD scheme both for 2-km and 500-m grid spacings. In the experiments starting at 0000 JST, the positional gap of the rainband with MYNN was smaller than that of the original experiments starting at 2100 JST owing to the short forecast time.
Similarly, in the experiments starting at 1800 JST, owing to the long forecast time, the positional gap in the MYNN experiments became larger, and the representation of the rainband deteriorated. Conversely, in the experiment with DD, the structure of the rainband was relatively well maintained even when the forecast time became longer. In addition, the horizontal distribution of the temperature in the lower layer and the updraft structure near the front for the experiments with different initial times was investigated. In the experiments starting at 0000 JST, the experiments with DD formed distinct vertical flows near the front earlier than the experiments with MYNN. In the experiments with both schemes, the fronts were maintained near the island and brought intense rain to the island.
In the experiments starting at 1800 JST, the experiments with DD formed distinct vertical flows along the front southeast of the island, and the fronts slowly shifted nearer to the island and were maintained. The experiments with MYNN also formed distinct vertical flows along the front southeast of the island; however, the front shifted northwest of the island and was not maintained near the island.
The difference between the experiments starting at 0000 JST and the experiments starting at 1800 JST was that the fronts were maintained near the island in the former experiments. Similar features were seen in the experiments starting at 2100 JST, as indicated in Section 4.2 (Figs. 9 -11). These results suggest that the DD scheme maintained the front near the island and brought intense rain. Meanwhile, the MYNN scheme did not maintain the front near the island, and the fronts shifted from southeast of the island to northwest. The impacts of the schemes were reproduced in the different initial time experiments.

Precipitation verification
In this section, we quantitatively evaluate the prediction skill of the precipitation distribution. For the verification method, the Fractions Skill Score (FSS) proposed by Duc et al. (2013) that extends the original score to both the spatial and temporal dimensions is used. The verification was conducted for the square domain E shown in Fig. 5. Figure 19 shows the FSS intensity-scale diagrams for KF5kmMY, CM2kmMY, CM2kmDD, CM500m MY, CM500mDD, and CM250mDD. The verification period was from 2300 JST October 15 to 0500 JST October 16. For a fair comparison, JMA's precipitation analysis data were up-(down-) scaled to model grids with grid spacings of 5 and 2 km (500 and 250 m), respectively. In this figure, the vertical axes display the spatial scales, which are the diameters of circular neighborhoods centered at each verification pixel. The bottom horizontal axes display the thresholds of the one-hour accumulated rainfall. Note that, while high thresholds consider only heavy rain, low thresholds take into account all events from light to heavy rain. The numbers on the top horizontal axes indicate the temporal scales in hours at each intensity threshold, e.g., 01 indicates no time lag, and 03 indicates a onehour time lag before and after the valid time.
In KF5kmMY (Fig. 19a), the FSS for the threshold of 1 mm h −1 was over 0.97 and that for the threshold of 10 mm h −1 was over 0.6. This indicates that KF 5kmMY was skillful in reproducing light and moderate rain (less than 10 mm h −1 ). When the threshold is over 30 mm h −1 , it is clear that the FSSs at the spatial scales of 10 and 15 km are much less than those at the spatial scales of 30 and 40 km. This results from the displacement of the intense rainband around Izu Oshima simulated by KF5kmMY. The same conclusion can be drawn from the FSSs calculated from the CM2kmMY experiment (Fig. 19b). The difference in the FSSs between CM5kmMY and CM2kmMY in Fig. 20a shows that these two experiments were comparable in predicting the precipitation, except for the slightly better skill of CM2kmMY at the threshold of 20 mm.
However, at the same resolution of 2 km, a significant difference can be observed with CM2kmDD (Fig. 19c), where the FSSs for the thresholds over 30 mm h −1 range from 0.20 to 0.60, which are much greater than the corresponding values observed with CM2kmMY. Furthermore, the positive values of the differences between the FSSs of CM2kmDD and CM2kmMY can be identified for thresholds ≥ 10 mm (see Fig. 20b), which indicates that CM2kmDD performed better than CM2kmMY. This improvement was achieved through the reduction in the positional lag of the rainband. This improvement was also confirmed at a resolution of 500 m when we compare the FSSs from CM500mMY (Fig. 19d) with those from CM500mDD (Fig. 19e) and their differences (Fig.  20c). Again, the reason for these differences was the reduction in the positional lag of the rainband when the DD scheme was used. Figure 19 also shows that a higher resolution tends to improve the rainfall forecast. Even though this improvement was nearly neutral in the MYNN experiments when we increased the resolution from 5 to 2 km, as depicted in Fig. 20a, a distinct improvement can be seen when we proceed from 2 km to 500 m in the DD experiment. Figure 20d shows that CM500m DD outperformed CM2kmDD for thresholds ≥ 10 mm. However, when the resolution was further increased from 500 to 250 m (CM250mDD in Fig. 19f ) with the same model orography, the FSSs slightly decreased (Fig. 20e).
To check the robustness of the verification results with respect to the verification domain size, the verification domain was zoomed into a square region centered on Izu Oshima with 180 km on each side. As shown in Supplements 4 and 5, all preceding conclusions still hold. Figure 21 shows the differences (a) between CM 500mDD (Fig. 8f ) and CM500mDD_small (Fig.  15e), (b) between CM500mDD and CM500m_nest (Fig. 15f ), and (c) between CM250mDD (Fig. 8g) and CM250mDD_K (Fig. 17d). Figure 21a, the score of CM500mDD was clearly higher than that of CM 500mDD_small for thresholds ≥ 10 mm. In Fig. 21b, the scores of CM500mDD for the thresholds of 20 -40 mm h −1 were slightly better than those of CM500m DD_nest. In Fig. 21c, the scores of CM250mDD_K were higher than those of CM250mDD for the intense rain area of ≥ 40 mm. These results indicate that the experiments with the large domain have better scores than the experiments with the small domain and that the experiment with the finer terrain improved the scores for intense rain.

Summary and concluding remarks
We performed ultra-high-resolution numerical weather prediction experiments using the JMA nonhydrostatic mesoscale model to examine whether the ultra-high-resolution model (500 to 250-m grid spacing) with a large domain (1600 km × 1100 km) can predict heavy rainfall more accurately. The Izu Oshima heavy rain event on October 16, 2013, was chosen. The JMA nonhydrostatic model optimized for the K computer was used to conduct these huge computations, and the computational efficiency was improved by approximately 10 -15 % on the K com-puter.
The impacts of (1) the grid spacing and (2) the PBL schemes [Mellor-Yamada-Nakanishi-Niino (MYNN) and Deardorff (DD)] were first investigated. The boundary layer scheme greatly influenced the position of the front and the associated rainband. In the experiments with DD, the front was simulated at a similar position to the observations, whereas in the experiments with MYNN, the front was shifted to the northwest of Izu Oshima. Convective clouds were generated with updrafts on the southeast side of the front, and downdrafts associated with rainfall occurred on the northwest side of the front. The experiments Vertical profiles of the model quantities (the potential temperature, water vapor, equivalent potential temperature, and wind speed) were examined on the southeastern (warm air) and northwestern (cold air) sides of the front. On the inflow side, vertical profiles of the model quantities of the warm and moist air simulated by the experiment with MYNN (CM500mMY) were slightly more upright than those simulated by the experiment with DD (CM500mDD) at the level from 300 m to 1.3 km. These differences are likely caused by the tendency of the MYNN scheme to yield larger vertical mixing in the boundary layer and the DD scheme to maintain the convective instability near the surface. On the cold air side, the heat flux at the sea surface was positive. The low-level potential temperature in CM500mDD was lower than that in CM500mMY; therefore, the temperature difference between the cold air from the Kanto Plain and the warm and moist air on the southeastern side of the front was larger than that in the experiments with MYNN. The vertical profiles of MYNN indicate that stronger vertical mixing occurred between the cold air in the lower level and the warm air aloft than occurred in experiments with DD. These differences also likely affected the intensity and the position of the front in each experiment.
The robustness of these findings was confirmed by experiments with different initial times. Similar responses to the two PBL schemes were obtained, with the positional gap of the simulated and observed rainbands becoming larger and the shape of the rainband To see the impacts of the model domain size on the rainfall prediction, experiments with the small model domain and nested simulations with different LBCs settings were examined. The experiment with a simple small domain exhibited a very different appearance of the rainband to that with the large domain. With the nesting procedure, the experiment with the small domain could reproduce a similar rainfall distribution to the large domain experiments if the cloud microphysical quantities were included in the LBC and the frequency of the boundary update was adequately high (every hour in our case). Note, however, that it is generally difficult to foresee the occurrence of a local heavy rainfall in real time NWP; therefore, such a very high-resolution simulation using nesting is only available after an event.
The influence of the terrain representation on the precipitation distribution was investigated. The 250-m grid experiment, with the finest terrain data generated by KTOPO, better reproduced the intense rain area on the island.
Performances of the quantitative precipitation forecast by the experiments with different horizontal grid spacings were verified using extended FSSs. In this case study, CM500mDD showed the best performance. CM500mDD with the large domain showed a slightly better performance than CM500mDD_ nest. The experiment with finer terrain representation (KTOPO) data improved the intense rain score in the verification.
In this case study, the 2-km grid experiment using DD simulated the position of the rainband better than what the experiment with MYNN did. The possible reason for this is that the DD scheme expressed the successive formation of new convection cells generated by the cold pool associated with the convectivescale updraft and downdraft at the front more distinctively than what the experiment with MYNN did. FSS verification for the 2-km grid-spacing experiments showed the advantage of DD for the prediction of very intense rains; however, the scores for weak to moderate rains (less than or equal to 10 mm h −1 ) were not improved.
In general, the direct application of the Deardorff scheme is only justified for very high resolutions, such as LES. At a grid spacing of 2 km, the MYNN scheme is usually applied (note that without boundary layer parameterization a large bias may occur in longterm forecasts). However the results here show that DD yielded better precipitation forecasts than MYNN at this grid spacing. An important research topic is therefore the development of an appropriate turbulence closure model for the "gray zone" (1 -2-km grid spacing). Several new planetary boundary schemes for the TI resolution have been proposed (e.g., Ito, J. et al. 2015;Kitamura 2016); however, most studies have been conducted for dry conditions, and tests of new TI schemes for heavy rain cases are future projects.
The overall results obtained in this study indicate that a very high-resolution model with a large domain has an advantage when simulating torrential heavy rain events; however, such a simulation requires huge computational resources. One solution to save computational resources is adaptive two-way nesting, which enhances the model resolution posteriorly by monitoring the prognostic valuables. Demonstrating the availability of such a strategy for heavy rainfall events is a future project.
Our experiments were conducted for a case study on a frontal rainband associated with a typhoon. Even though the results were robust for different initial times, the situation is likely not the same in other cases. Recently, the authors investigated another heavy rainfall event in Hiroshima on August 19 -20, 2014, whose horizontal scale was much smaller than the Izu Oshima heavy rain event. In that case, different factors determined the rainfall position, and the detailed results will be shown in another paper.
Development of a high-performance computing debris flow model is underway, and the direct prediction of disasters by a hydrological model using forecasts of the ultra-high-resolution NWP is an important research subject. Such studies will be performed in the FLAGSHIP2020 project using the Post K computer.

Supplements
Supplement 1 shows the sensible heat flux at the model surface of (a) CM500mMY and (b) CM500m DD at 0100 JST on October 16. Supplement 2 is the same as Figs. 13a -d except for CM500mDD_D'. Panel (e) shows the sampling area for CM500mDD_ D'. Supplement 3 is the same as Figs. 15a -d except that the model domain size is a 400-km 2 region. Supplement 4 is the same as Fig. 19 except for the small verification domain size. Supplement 5 is the same as Fig. 20 except for the small verification domain size.

Acknowledgments
A part of this work was supported by the Ministry of Education, Culture, Sports, Science and Technology as the Field 3, the Strategic Programs for Innovative Research (SPIRE), the FLAGSHIP2020 project (Advancement of meteorological and global environmental predictions utilizing observational "Big Data"), and the JSPS Grant-in-Aid for Scientific Research 'Study of optimum perturbation methods for ensemble data assimilation' (16H04054). The authors thank Dr. Hiromu Seko of the Meteorological Research Institute for his continuous support and encouragement.
Computational results were obtained using the K computer at the RIKEN Advanced Institute for Computational Science (project ID: hp140220, hp150214, and hp160229). A part of input and output data were processed with UV1000 at the Japan Agency for Marine-Earth Science and Technology. Detailed and constructive comments on the manuscript by two anonymous reviewers and the editor, Dr. Masayuki Kawashima of the Hokkaido University, significantly improved the quality of this paper. We give thanks to Dr. Akira Noda of Japan Agency for Marine-Earth Science and Technology for providing useful comments.