Flood Simulations in Mid-Latitude Agricultural Land Using Regional Current and Future Extreme Weathers

: Recent extreme weather events like the August 2016 ﬂood disaster have signiﬁcantly a ﬀ ected farmland in mid-latitude regions like the Tokachi River (TR) watershed, the most productive farmland in Japan. The August 2016 ﬂood disaster was caused by multiple typhoons that occurred in the span of two weeks and dealt catastrophic damage to agricultural land. This disaster was the focus of our ﬂood model simulations. For the hydrological model input, the rainfall data with 0.04 ◦ grid space and an hourly interval were provided by a regional climate model (RCM) during the period of multiple typhoon occurrences. The high-resolution data can take account of the geographic e ﬀ ects, hardly reproduced by ordinary RCMs. The rainfall data drove a conceptual, distributed rainfall–runo ﬀ model, embedded in the integrated ﬂood analysis system. The rainfall–runo ﬀ model provided discharges along rivers over the TR watershed. The RCM also provided future rainfall data with pseudo-global warming climate, assuming that the August 2016 disaster could reoccur again in the late 21st century. The future rainfall data were used to conduct a future ﬂood simulation. With bias corrections, current and future ﬂood simulations showed the potential inundated areas along riverbanks based on ﬂood risk levels. The crop ﬁeld-based agricultural losses in both simulations were estimated. The future cost may be two to three times higher as indicated by slightly higher simulated future discharge peaks in tributaries.


Introduction
Future climate changes may significantly affect the mid-latitude region worldwide according to reports by the Intergovernmental Panel on Climate Change (IPCC) [1,2]. The typhoon-affected areas in the North Pacific have been potentially extended to the mid-latitude region. However, past studies related to typhoon-impact on flood disasters in the mid-latitude region are relatively few as flood events are rare in the region. For example, academic reports and articles in English on typhoon-induced flood disasters have seldom been published in the northern East Asian regions such as Northeast China and the northern Korean Peninsula (e.g., Myeong et al. [3]). We focused on a mid-latitude area in the North Pacific, Hokkaido, Japan's northern island, ranging between 41. 35-45.56 • N and 139. .89 • E. Hokkaido has recently experienced several severe flood events due to typhoon-included heavy rainfall, although typhoons were rare between 1961 and 2016, according to the Japan Meteorological Agency  (1) The data were recorded when the typhoons landed or most closely approached Hokkaido, and are referred to as the Best Track dataset in the RSMC Tokyo [4]. (2) The data were observed at Nukabira-Gensenkyo Station in the eastern portion of Hokkaido [8].

Integrated Flood Analysis System (IFAS)
We performed several numerical flood examinations in the main and tributary rivers over the TR watershed. We employed a comprehensive, GUI-operating runoff simulator, the IFAS, combined with satellite weather data. The IFAS was developed by the International Centre for Water Hazard and Risk Management (ICHARM) [18,19], and has been applied to historical flood events all over the Kitano et al. [9] showed that some typhoons after 2011 tracked across Hokkaido on different courses when compared with past typhoons due to a continuous stay of the same atmospheric pressure distributions. This phenomenon potentially prompted the occurrences of heavy rainfall events. Recent studies related to the statistical analyses of future climate projections in Hokkaido [10,11] reported that future heavy rainfall events and related-flood risks will be increased as per ensemble simulations based on the database for policy decision making for future climate change (d4PDF), which involves past climate data (3000-year numerical experiments) and future 4 • C increase-based climate data (5400-year numerical experiments) [12]. Therefore, future flood disasters over Hokkaido may be more severe. On the other hand, climate studies regarding specific typhoons have also been conducted. For example, a regional climate model (RCM) was recently run to analyze current and future climates for specific typhoon tracks that caused the severe disaster in Hokkaido in August 2016 [13].
Unlike the statistical analyses of flood disasters conducted using the d4PDF database, we focused on the worst flood disaster in the region, caused by the multiple typhoons during late August 2016 (hereafter, the August 2016 disaster), using regional reproduction of current climate and future climate data. We analyzed the flood-damaged farmland, spread in the southwest portion of the Tokachi River (TR) watershed (42.54-43.64 • N, 142.68-144.03 • E), and located in the southeast region of Hokkaido ( Figure 1a). Heavy rainfall events occurred during August 16-24 when the preceding three typhoons struck (e.g., >400 mm at the Naitai Station), according to the water resource database provided by the Ministry of Land, Infrastructure, Transport, and Tourism of Japan (MLIT Japan) [8]. The fourth typhoon, Typhoon Lionrock, on August 29-31 passed the south of Hokkaido and brought heavy rainfall (e.g., >200 mm at the Naitai Station) around the Daisetsu Mountains (central Hokkaido) and over the westside of the Hidaka Mountains (southeast Hokkaido). The meteorological data suggested that the raincloud hit the Hidaka Mountains from the west, and then the upward movement of air brought approximately 500 mm of precipitation for 48 h at the Satsunai River Dam Station near the Hidaka Mountains. In particular, severe floods occurred in several tributaries of the southwest portion of the TR watershed, thus damaging agricultural areas [14]. As the agricultural areas are the most productive farmland in Hokkaido, the economic losses were excessive [7].
In other research related to the August 2016 disaster, Kimura et al. [15] showed that with a hydrological model and observed rainfall data, the effect of the preceding three typhoons on flood disasters was crucial as the water capacity in soil in the west side of the TR watershed was exceeded before the fourth typhoon approached, which then made the flood disaster worse by producing massive rainfalls. However, this previous study was limited to understanding the future trend because there is no future rainfall data based on climate change scenarios as an extended investigation of the August 2016 disaster. In the recent domestic situation for agricultural losses, major-typhoon impacts on the losses of agriculture, forestry, and fisheries in Japan have been more severe than a decade ago. The annual losses were roughly estimated as over a billion USD from 2016 to 2018 [16]. Therefore, to maintain sustainable agricultural activities in the TR watershed, it is necessary to evaluate future flood risks over the farmland. We will show the difference in agricultural losses between the actual-disaster and possible-flood damages under the impact of climate change by using an RCM.
This study focused on evaluating the crop-field-based agricultural damage loss over the inundated areas in the southern portion of the TR watershed during the August 2016 disaster using current and future extreme weather data provided by the RCM. The new trials were to evaluate the August 2016 disaster with the reproduction current and future climates and estimate the agricultural damage costs in the specific future climate as a first trial. This study is important to reveal how much of the agricultural damage may be changed quantitatively in a specific area as an illustrative example. Note that we employed the flood model, appropriately calibrated in our previous study using on-site observations [15].

Study Site
The TR watershed, located in the eastern portion of Hokkaido, has an area of approximately 9010 km 2 , stretching along approximately 156 km of the Tokachi River through low-plains to the North Pacific with a maximum altitude of 2077 m (Figure 1b). Most of the low-plain areas are utilized as farmland. The tributaries of the Tokachi River are still vulnerable to severe flooding because of poor countermeasures, although a few dams and weak embankments for flood control exist along the tributaries [17]. Regarding the hydraulic features of the watershed, heavy rainfall events usually occur in the high mountains, which stand at the west, north, and east boundaries of the TR watershed; the rainwater gathers downstream through the river network on the right and left sides of the watershed and flows to the outlet in the North Pacific. A time lag of several hours exists before heavy rainfall in the upstream area reaches the downstream area and the outlet.

Integrated Flood Analysis System (IFAS)
We performed several numerical flood examinations in the main and tributary rivers over the TR watershed. We employed a comprehensive, GUI-operating runoff simulator, the IFAS, combined with satellite weather data. The IFAS was developed by the International Centre for Water Hazard and Risk Management (ICHARM) [18,19], and has been applied to historical flood events all over the world [20][21][22]. A conceptual, distributed rainfall-runoff analysis engine, Public Works Research Institute-distributed hydrological model (PWRI-DHM) [23], is embedded in the IFAS as a module. Several coefficients in the PWRI-DHM are efficiently set up based on soil-type and land use features. Detailed explanations are shown in Appendix A. For the IFAS computational conditions in the TR watershed, the number of cells was 170 × 137 in the horizontal direction over the entire watershed and each cell size was 0.008 • × 0.008 • . The river channel network over the watershed was determined using the altitude information.

Data Inputs
The IFAS required the elevation, land use, and soil-geology data to be used as hydrological information inputs. The elevation data were obtained from the digital elevation model that used altitude data from the shuttle radar topography mission with a three arc-second resolution (approximately 90 m) in the hydrological data and maps based on shuttle elevation derivatives at multiple scales [24]. The land use data with a 1 km mesh resolution collected through the U.S. Geology Survey were downloaded from the Global Land Cover Characterization Landuse website [25]. The data for the geology and soil type were obtained from the global distribution data for soil water holding capacity with 1 • resolution between 0 to 0.3 m in the surface of the soil. The global distribution data were compiled as part of the United Nations Environment Program. The rainfall data for the model inputs were provided by the RCM and then distributed into the separate watershed divisions using the ordinary kriging, which is a geostatistical interpolation technique using a least-squares regression method among neighborhood observations.

Rainfall Data and Bias Corrections
A high-resolution, three-dimensional atmosphere-ocean coupled regional model, the Cloud Resolving Storm Simulator version 3.4 (CReSS) [26,27] for the atmospheric part and the Non-Hydrostatic Ocean model for the Earth Simulator (NHOES) [28,29] for the oceanic part (hereafter CReSS-NHOES), provided the rainfall datasets during multiple typhoons under the conditions of controlled current (CNTL) climate as a present simulation and pseudo-global warming (PGW) climate as a future simulation [13]. The CReSS-NHOES simulated the horizontal grid spacing, ranging from 132 • E to 155 • E and 25 • N to 50 • N, discretized with a 0.04 • × 0.04 • grid, and hourly temporal resolution [13]. First, the CNTL simulation was conducted. Initial and lateral boundary conditions for the atmosphere and ocean models were provided by the Japan Meteorological Agency 55-year Reanalysis dataset [30] and the Japan Coastal Ocean Predictability Experiment reanalysis product (JCOPE2) [31], respectively.
Second, we created future climate data (hereafter, PGW simulation) considering the climate changes between current and future climate conditions. The climate changes were made using the atmospheric general circulation model of the Meteorological Research Institute (MRI-AGCM) version 3.2 run with a horizontal resolution of 20 km [32]. A multi-model ensemble mean sea surface temperature (SST) pattern projected by Coupled Model Inter-comparison Project phase 5 models under the representative concentration pathway (RCP) 8.5 scenario with a 4 • C increment [33] was employed as the computational condition of future climate. Subsequently, the changes in the August monthly mean between the periods of 1979−2003 and 2075−2099 were added to the initial and boundary conditions of the CNTL on the variables of the SST, atmospheric temperature, and water vapor.
The regional rainfall datasets in the TR watershed during the multiple typhoons in August 2016 were obtained from the rainfall datasets in the CNTL and PGW simulations. The accumulated quantities of both rainfall datasets are shown with the observed precipitation data in Figure 2. A comparison between the two datasets shows that the CNTL simulation and observed data were similar, although the CNTL simulation was slightly underestimated. Figure 3 shows the maximum intensities per hour among the CNTL, PGW, and observed data. The CNTL simulation was much higher than the observed data. In addition, the accumulation quantity for the PGW simulation was lower than that of the CNTL simulation, but the PGW simulation had a maximum intensity.  Due to the differences between the CNTL simulation and observed data, we modified the CNTL simulation using two bias correction methods: a simple linear method and a statistical method [34] based on hourly spatial rainfall data of approximately 50 gauge stations covering the TR watershed ( Figure 1b). The observed rainfall data were obtained from MLIT Japan and JMA [8,35]. The linear method adjusts the CNTL simulation to the observed data using a linear relationship between both sets of precipitation data. Both sets of hourly rainfall data during the target period were arranged in ascending order of quantity, before the regression line was generated. In contrast, the statistical method adjusts the CNTL simulation using cumulative distribution functions (CDF) of both datasets (CNTL simulation and observed data). However, the data correction methods cannot modify the timings at which the highest amount of rainfall occurs because they only adjust the total amount of precipitation between the observed and projected data during the target period. Furthermore, the PGW simulation was modified using a combination of the linear and statistical bias correction methods. Detailed explanations are shown in Appendix B. Figure 2. Accumulated distributions of the precipitation during the August 2016 disaster. "cntl", "pgw", and "observation" indicate reproductions of the current climate, future climate, and real world observed data at rain-gauge stations, respectively.

Simulation Preparations
We employed the PWRI-DHM for the flood simulation. The model was calibrated in our previous study using appropriate parameters, along with the spatial and temporal resolutions [19]. Our previous work presents the flood model setup in detail [15]. Flood control dams were implemented in the flood model. The flood risk discharges were estimated using the rating curve between the discharge and water stages, derived from past observed data from each gauge station using available risk level information (i.e., a bank-full water level) on the water stage [8]. The two major dams, Tokachi Dam and Satsunai River Dam, used for flood control (Table 2) were set up for the model. However, their flood control roles only involved retaining an amount of water based on the dam capacities as data for flexible flood controls such as adaptive control (e.g., pre-release discharge) and abnormal-flood-prevention operations could not be estimated using the climate simulations, unlike the dam flood controls that were conducted in reality. Note that the model did not consider flow loss from rivers due to broken dykes as the flooding areas were within a one grid scale, and an uncertainty analysis from observation designs such as the number of rain-gauge stations and the appropriate locations of these stations was beyond the scope of this study.

Indicators for Model Validation
To evaluate how accurately the flood model simulation matched the observed data, we introduced the Nash-Sutcliffe model efficiency coefficient (NS) [36] for a whole wave-shape reproducibility and the mean absolute error (MAE) for wave peak reproducibility. The NS is the following equation.
where V O is the observed data; V M is the simulated data; N is the number of points; and V AVG is the average value of the observed data. If the value of NS is over 0.7, the simulation reproducibility can be accepted as a quantitative evaluation. The MAE is also an effective indicator to measure the error, which is given by Note that hourly collected discharge data were obtained from MLIT Japan [8].
All of the procedures based on data flow are depicted in Figure 4 in the integrated methodology of this study. where O V is the observed data; M V is the simulated data; N is the number of points; and AVG V is the average value of the observed data. If the value of NS is over 0.7, the simulation reproducibility can be accepted as a quantitative evaluation. The MAE is also an effective indicator to measure the error, which is given by Note that hourly collected discharge data were obtained from MLIT Japan [13]. All of the procedures based on data flow are depicted in Figure 4 in the integrated methodology of this study.

Results
We validated the reproduction current flood simulation by the hydrological model applied to the TR watershed during the time period between August 16 and 5 September 2016, which indicated the multiple typhoons. The validation was conducted to compare the current flood simulation with the observed discharge data obtained from several gauge stations (Figure 1b), where agricultural damage was clustered. We tested two bias correction methods before conducting the future flood simulation during the same period with the PGW simulation data. Furthermore, we estimated the crop-field-based agricultural damage losses for the August 2016 disaster. Note that severe overflows occurred in several tributaries in the west and southwest portions of the TR watershed such as the Pekerebetsu River and Memuro River. However, observed datasets of the discharge were limited, so we only considered the following tributaries with available datasets: Bisei River, Tottabetsu River, and Satsunai River (Table 3).

Results
We validated the reproduction current flood simulation by the hydrological model applied to the TR watershed during the time period between August 16 and 5 September 2016, which indicated the multiple typhoons. The validation was conducted to compare the current flood simulation with the observed discharge data obtained from several gauge stations (Figure 1b), where agricultural damage was clustered. We tested two bias correction methods before conducting the future flood simulation during the same period with the PGW simulation data. Furthermore, we estimated the crop-field-based agricultural damage losses for the August 2016 disaster. Note that severe overflows occurred in several tributaries in the west and southwest portions of the TR watershed such as the Pekerebetsu River and Memuro River. However, observed datasets of the discharge were limited, so we only considered the following tributaries with available datasets: Bisei River, Tottabetsu River, and Satsunai River (Table 3). (1) rMAE is defined by MAE/observed discharge peak × 100 (%).

Hydrological Model Validation
The validation of the hydrological model with and without a bias correction was conducted. Using the raw data over the TR watershed from the CNTL simulation, the simulated time-series discharges are shown ( Figure 5) at eight gauge stations, where Moiwa Station (#1 in Table 3), close to the end of the main stream, was added as a reference to check the entire discharge over the TR watershed. The other stations are located along tributaries vulnerable to flood disasters. In particular, discharge peaks in late August were underestimated, although most wave shapes were similar to the observed ones as indicated by the NS values (Table 3). However, the discharge peaks were unable to be appropriately captured. The accuracy of the peak discharge is more important than that of the wave shape of the discharge because our study aimed to calculate the overflow volume when the simulated discharge was greater than the flood risk discharge. The temporal discharge at Moiwa Station was overestimated during mid-August. This is because the Tokachi Dam, located on the upper-side of the main stream, was not effective at reducing the flood risk in the downstream areas as adaptive flood control was not included in the flood model.
The simulated result provided by the linear relationship between the CNTL simulation and observed data in the ascending coder was named, "Linear BC." The simulated discharge at Moiwa Station was similar to the observed one with a NS of 0.70 (Table 3), although the maximum peak was not captured appropriately during the approach of Typhoon Lionrock on August 29-31 (Figure 6a) due to less estimated precipitation than the observation over the whole TR watershed (Figure 2). The temporal discharges at the other stations in the tributaries were good with NSs of approximately 0.70 and rMAEs between 2.7 and 5.2 in comparison with the observed data, although the discharge peaks were not reproduced well (Figure 5b-h). Note that the maximum peaks at some stations such as Totta Bridge and Nakajima Bridge Stations seem to be missing; however, the observed data at the other stations along tributaries, that did not have any missing data, showed that the peaks dropped suddenly after the maximum discharges were recorded. Thus, the missing data might not have a significant impact on the validation of the model. The simulated result provided by the linear relationship between the CNTL simulation and observed data in the ascending coder was named, "Linear BC." The simulated discharge at Moiwa Station was similar to the observed one with a NS of 0.70 (Table 3), although the maximum peak was not captured appropriately during the approach of Typhoon Lionrock on August 29-31 ( Figure 6a) due to less estimated precipitation than the observation over the whole TR watershed (Figure 2). The temporal discharges at the other stations in the tributaries were good with NSs of approximately 0.70 and rMAEs between 2.7 and 5.2 in comparison with the observed data, although the discharge peaks were not reproduced well (Figure 5b-h). Note that the maximum peaks at some stations such as Totta Bridge and Nakajima Bridge Stations seem to be missing; however, the observed data at the other stations along tributaries, that did not have any missing data, showed that the peaks dropped suddenly after the maximum discharges were recorded. Thus, the missing data might not have a significant impact on the validation of the model. Furthermore, we tested "Statistical BC", the other bias correction method. As seen in Figures 7bh, the stations along the tributaries gained better discharge peaks than the peaks obtained through Linear BC and relatively higher NSs values than those of Linear BC, except for the Moiwa and Dainiohkawa Bridge Stations. The magnitudes of rMAEs were relatively lower than, or approximately equivalent to, those of the other two simulations, which implied that the discharge Furthermore, we tested "Statistical BC", the other bias correction method. As seen in Figure 7b-h, the stations along the tributaries gained better discharge peaks than the peaks obtained through Linear BC and relatively higher NSs values than those of Linear BC, except for the Moiwa and Dainiohkawa Bridge Stations. The magnitudes of rMAEs were relatively lower than, or approximately equivalent to, those of the other two simulations, which implied that the discharge peaks were quantitatively more accurate. Overall, considering the discharge wave-shape and peak, the flood simulation with Statistical BC showed good agreement with the observed data. Therefore, from the above three simulation results on the CNTL simulation, we can see that the results from Statistical BC was a better representation of the simulated discharge and could appropriately capture the peaks.

Future Flood Situation under GLOBAL WARNING
The future flood simulation conducted using the PGW simulation data and a combination of Statistical BC and Linear BC was compared with the current flood simulation conducted with Statistical BC. Figure 8 indicates the discharge differences between the CNTL and PGW simulations at the eight stations. Simulated future discharges at seven stations along the tributaries were 2-3% larger than those of the current flood simulation (Figure 8b-h). The future flood discharge at Moiwa Station was different from the current discharge with an approximately 7% reduction (Figure 8a).

Future Flood Situation under GLOBAL WARNING
The future flood simulation conducted using the PGW simulation data and a combination of Statistical BC and Linear BC was compared with the current flood simulation conducted with Statistical BC. Figure 8 indicates the discharge differences between the CNTL and PGW simulations at the eight stations. Simulated future discharges at seven stations along the tributaries were 2-3% larger than those of the current flood simulation (Figure 8b-h). The future flood discharge at Moiwa Station was different from the current discharge with an approximately 7% reduction (Figure 8a).

Estimation of Agricultural Damage Loss
We estimated the agricultural damage losses based on the overflow volumes from the flood simulations (Statistical BC) and the production costs for the agricultural areas with several assumptions such as overflow occurrences at gauge stations, uniform blocks over crop fields, and no erosion processes due to limited datasets and little information on the processes. For example, there were numerous locations that experienced overflow events in reality. However, the observed datasets were not available in those overflow locations. Instead, we considered the possible damage losses only at several gauge stations that did not actually experience overflows, but had the observed discharge data available. This is because the purpose of our study was to understand the difference in agricultural damage caused by the impacts of current and future climates.
According to the annual report of the statistical data at the Hokkaido Agriculture, Forestry, and Fisheries [37] and our cost calculations on the major products for latest five years (Table 4), the average profit of regional agricultural products such as potatoes, beans, and sugar beet was estimated to be 0.72 M-USD per km 2 . We only focused on the agricultural product damages because it is difficult to estimate the damage costs for facilities and machines on farmland due to the need to consider the individual damaged situations. We assumed that the production profit per unit area was equivalent to the cost of damage caused by the flood disaster. We also assumed that the overflow volume was obtained at each gauge station only when the flood simulation discharge (qr, m 3 /s) was greater than the flood risk discharge (qc, m 3 /s). The overflow volume was calculated by ∑ ( , − , ) × , where i is a numbering of time step, and DTi is the duration (3600 s) when the overflow occurred (i.e., qr,i -qc,i > 0). We then calculated the inundated area by dividing the overflow volume by the height (approximately 0.4 m) from the ground on a typical crop field to the road along the field in which the agricultural products are usually destroyed. This calculation did not consider multiple overflow

Estimation of Agricultural Damage Loss
We estimated the agricultural damage losses based on the overflow volumes from the flood simulations (Statistical BC) and the production costs for the agricultural areas with several assumptions such as overflow occurrences at gauge stations, uniform blocks over crop fields, and no erosion processes due to limited datasets and little information on the processes. For example, there were numerous locations that experienced overflow events in reality. However, the observed datasets were not available in those overflow locations. Instead, we considered the possible damage losses only at several gauge stations that did not actually experience overflows, but had the observed discharge data available. This is because the purpose of our study was to understand the difference in agricultural damage caused by the impacts of current and future climates.
According to the annual report of the statistical data at the Hokkaido Agriculture, Forestry, and Fisheries [37] and our cost calculations on the major products for latest five years (Table 4), the average profit of regional agricultural products such as potatoes, beans, and sugar beet was estimated to be 0.72 M-USD per km 2 . We only focused on the agricultural product damages because it is difficult to estimate the damage costs for facilities and machines on farmland due to the need to consider the individual damaged situations. We assumed that the production profit per unit area was equivalent to the cost of damage caused by the flood disaster. We also assumed that the overflow volume was obtained at each gauge station only when the flood simulation discharge (q r , m 3 /s) was greater than the flood risk discharge (q c , m 3 /s). The overflow volume was calculated by numbering of time step, and DT i is the duration (3600 s) when the overflow occurred (i.e., q r,i −q c,i > 0). We then calculated the inundated area by dividing the overflow volume by the height (approximately 0.4 m) from the ground on a typical crop field to the road along the field in which the agricultural products are usually destroyed. This calculation did not consider multiple overflow events (i.e., a single overflow event at each gauge station was assumed). Note that this kind of simple calculation is often used to gain the first cut estimation of flood damage costs based on land-use classification in Japan [38] with the damage rate being approximately 1.0. Average gross profit for the major products (JPY/10 3 m 2 ) 76,724 Average gross revenue for the major products (M-USD/km 2 ) 0.72 A typical agricultural cost in the TR watershed (0.72 M-USD/km 2 ) was calculated based on the agricultural information from the website of the Obihiro office [37]; the exchange rate of JPY to USD was assumed to be 0.0094; M stands for million.
The estimated damage losses between the current and future flood simulations are shown in Table 5. Except for Moiwa Station, the future damage costs were two to three times greater than the current damage costs at three stations: Nakajima Bridge, Rrutanjyoryu, and Minamisatsunai. The damage costs at Dainiohkawa Bridge and Nantai Bridge in the Satsunai River tributary are currently minimal, but the future flood simulation showed that damages of 3.9 and 0.4 M-USD are expected to occur at the respective stations. This indicates that flooding is expected to occur here in the future as a result of potential heavy rainfall events in the east side of the Hidaka Mountains.

Discussion
Our flood simulation trails were conducted to obtain a better representation of the temporal discharge during the August 2016 disaster (Figures 5-7). With the raw rainfall data from the CNTL simulation at the spatial resolution of 0.04 • × 0.04 • , the simulated discharges at several stations were quite similar to the observed ones based on the quantitative evaluation of the wave shape (i.e., NSs > 0.7 for the five stations in Table 3). This result possibly came from the local geographic effects being included in the precipitation. Hoshino et al. [10] suggested that the regional downscaling model with a high-spatial resolution (5 km, approximately equivalent to 0.04 • ) provides highly accurate precipitation estimates close to the observed extreme rainfall events. Other studies [39,40] also pointed out that spatially high-resolutions (2-4 km) in dynamic downscaling can be required to resolve the precipitation over complicated topography. Therefore, we can conclude that the local precipitation cut directly from the CNTL simulation was acceptable for the use of our flood simulation as general climate simulations for past or current events cannot be perfectly reproduced.
However, our study needed to calculate an overflow volume when the simulated discharge was beyond the flood risk discharge at each station. Thus, accurately simulating a discharge peak was crucial. We introduced the bias correction methods to improve the discharge peak accuracy. The flood simulation with the Statistical BC provided a better representation of the temporal discharge during the August 2016 disaster (Figure 7). This bias correction method modifies the rainfall distribution of the CNTL simulation into a more appropriate distribution based on the observed rainfall using the same CDF of rainfall intensities between the CNTL simulation and observation. Consequently, it was necessary to introduce Statistical BC to the modification of precipitation.
According to our previous study [15], the soil had exhausted its water capacity in the southwest portion of the TR watershed before the fourth typhoon and, as a result, the severe flood disaster occurred during the fourth typhoon's approach toward Hokkaido. The agricultural damage areas at the Nakajima Bridge, Rrutanjyoryu, and Minamisatsunai Stations calculated in the current flood simulation were 3-10 km 2 along the tributaries where overflow events occurred, but not at Moiwa Station. These estimated areas are reasonable as the actual damage area in the town close to Bisei Bridge Station was approximately 9.1 km 2 [14], although the damage was unrelated to agricultural products. The damage cost was 2.2-7.8 M-USD, approximately equivalent to the actual damage cost (2.8 M-USD) in the farmland where there are some tributaries in the southwestern TR watershed [41]. The damage cost was approximately 3-4% of the total damage cost in Hokkaido [7]. These data for local areas can be helpful to farmland stakeholders as the farmland damage costs have not usually been made public. The future damage costs to the farmland were two to three times larger than the current ones at the three stations. Moreover, similar damage costs occurred at two other stations: Dainiohkawa Bridge and Nantai Bridge. These increments were caused by the higher rainfall intensities of the PGW simulation around the southwestern TR watershed (Figure 3).
In addition, the future time-series discharge at Moiwa Station, close to the end of the main stream (i.e., the outlet of the river), was significantly different from the current discharge ( Figure 8a). The large difference can potentially be explained by the change in the temporal and spatial distributions of precipitation. Figure 9a shows the potential typhoon trajectories between the raw datasets of the CNTL and PGW simulations by plotting the lowest air pressure at each hourly time step around Hokkaido. In particular, the biggest typhoon, Typhoon Lionrock, had slightly different trajectories between both datasets from August 28 to 31. In fact, the time series precipitation averaged over the upper area from Moiwa Station showed that the rainfall of the PGW simulation was lower than the CNTL simulation with reduced values of approximately 20% over the entire period and 3% from August 28 to 31 (Figure 9b). For the change in the spatial distribution, the future precipitation was weaker in accumulated quantity in Figure 2 and had a higher intensity only around the southwestern TR watershed (Figure 3). between both datasets from August 28 to 31. In fact, the time series precipitation averaged over the upper area from Moiwa Station showed that the rainfall of the PGW simulation was lower than the CNTL simulation with reduced values of approximately 20% over the entire period and 3% from August 28 to 31 (Figure 9b). For the change in the spatial distribution, the future precipitation was weaker in accumulated quantity in Figure 2 and had a higher intensity only around the southwestern TR watershed (Figure 3).  Worldwide discussions about the impacts of future climate change on agriculture damage in nations and regions have been actively undertaken during the last two decades. For example, the damage estimations of agricultural products caused by climate change of the 2 • C-up target in the year 2100 were conducted using the environmental damage-adaptation model including economic factors such as labor and capital and damage states. The model provided the agricultural damage index, which shows the negative effect on grain and crop products at 2050 in Japan [42] in comparison with the 2004 index. The KAKUSHIN program reported that using a climate-ocean model incorporating the dynamics of ecosystem, chemical, and other processes, the crop yield (mainly maize) reduction rates estimated in the RCP 8.5 scenario for USA, Brazil, and China were increased during the period of 2040-2070 due to climate change [43]. The above studies suggest that the near-future climate could negatively affect agricultural products in the relatively large-scale impacts on agricultural damage. However, our study focused on the impact for a small-scale watershed. Our small-scale results had a similar trend of the negative impacts of climate change on large-scale agricultural damages. In the relatively small-scale watershed (21,734 km 2 ) in Turkey, the agricultural productivity change between the 2000s and 2070s was indirectly estimated based on the shift of the plant allocations on farmland using the atmospheric general circulation model's outputs and risk aversion classifications [44]. However, their results did not show specific agricultural losses by climate change. Agricultural damage over crop fields in a small-scale watershed have likely not been estimated. Therefore, our results of the agricultural damage by the specific flood event can be useful, at least as a first-cut estimation under future climate change impacts.

Conclusions
We estimated the possible agricultural damage losses in the southwest portion of the TR watershed, where a severe flood disaster occurred during late August 2016, by using the reproduction-current and future climate data. First, the reproduction flood simulations were tested using the CNTL simulation with and without bias correction methods (Linear BC and Statistical BC). The current flood simulation with Statistical BC was a better representation of the time-series discharge, considering the reproducibility of discharge peaks, with NS values in the range of 0.4-0.9 and rMAE values between 2.4-5.2% at seven stations along the tributaries. Second, the future flood simulation was performed with a combination of Linear BC and Statistical BC. The future discharges were slightly larger than the current ones. Finally, both flood simulations in the CNTL and PGW simulations provided possible overflow volumes as an illustrative example due to several assumptions. We then estimated the agricultural damage costs simply calculated by the overflow volumes, which suggested that future damage costs might be two to three times greater, possibly due to the higher intensities of precipitation in the southwestern TR watershed in the future climate.
In our future work, to reduce the uncertainties for future climates, the PGW simulations with possibly different scenarios and a higher-resolution model will be necessary as an extension of the single SST-conditional simulation. Furthermore, as the agricultural damage areas were simply calculated with several assumptions such as no erosion processes, uniform blocks over the crop field, and simplified overflow volumes due to limited datasets and little information of those, the next step requires taking into account the specific factors that may cause flooding over crop fields. Acknowledgments: Numerical simulations of the climate models were performed using the supercomputing system at Nagoya University. We greatly appreciate Kazuhisa Tsuboki (Nagoya University, Japan) for his helpful support.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A
As a flow simulation engine, the PWRI-DHM uses the tank model, which consists of surface, groundwater, and river tanks. The tank model concept is depicted in Figure A1. Surface flow (q sf ) and subsurface flow (q ri ) are determined by surface water level (sWL) and infiltration (q o ) in the surface tank; in the groundwater tank, unconfined groundwater flow (q g1 ) next to the groundwater tank and confined groundwater flow (q g2 ) to the river tank are estimated by ground water level (gWL) and input flow (q gin ) from the surface tank; river flow (q r ) is calculated from the difference between the water level (rWL) in the river tank and input flow (q rin ) from the surface and groundwater tanks. Further information is available on the ICHARM website [19]. tank; in the groundwater tank, unconfined groundwater flow (qg1) next to the groundwater tank and confined groundwater flow (qg2) to the river tank are estimated by ground water level (gWL) and input flow (qgin) from the surface tank; river flow (qr) is calculated from the difference between the water level (rWL) in the river tank and input flow (qrin) from the surface and groundwater tanks. Further information is available on the ICHARM website [19]. Figure A1. Schematic of the three-tank model in the public works research institute-distributed hydrological model.

Appendix B
As local precipitation from climate datasets needs to be adjusted to the observed data, we utilized the linear and statistical bias correction methods between the CNTL simulation and observed data. Figure A2a shows the concept of the linear bias correction. Both precipitation datasets were arranged individually in the ascending order of quantities. Then, the regression line between both datasets was generated. Figure A2b describes how to use the statistical bias correction, based on observed CDF, to adjust the CNTL simulation. The principle is that the CDF of the CNTL simulation should be the same as that of the observed data to adjust the total amount of precipitation. Although this bias correction has usually been used for long-term daily precipitation, we assumed that hourly data within a short-term period (approximately two weeks) could be applied as limited data were available. Note that in general, it is difficult to use the statistical bias correction to adjust the timing at which peak precipitation intensities are observed. Furthermore, the bias correction of the PGW simulation was conducted in two steps. We first conducted the statistical bias correction to which the parameters in the CNLT bias correction were simply applied and then further modified the corrected data using the linear relationship between the CNTL and PGW simulations to consider the trend of future climate change based on a similar idea by Pierce et al. [45] (Figure A2c). Figure A3 shows the accumulated amounts of precipitation among the observed and two bias-corrected precipitations for the CNTL simulation during the two weeks. The Linear BC distribution of the accumulated precipitation was approximately equivalent in amount to the observed distribution. However, the Statistical BC distribution was 10% higher than the observation. This suggests that the observed distribution could be underestimated because the kriging interpolation was performed with sparser data points than the grid spacing of the CNTL simulation, especially in the east portion of the TR watershed. Figure A1. Schematic of the three-tank model in the public works research institute-distributed hydrological model.

Appendix B
As local precipitation from climate datasets needs to be adjusted to the observed data, we utilized the linear and statistical bias correction methods between the CNTL simulation and observed data. Figure A2a shows the concept of the linear bias correction. Both precipitation datasets were arranged individually in the ascending order of quantities. Then, the regression line between both datasets was generated. Figure A2b describes how to use the statistical bias correction, based on observed CDF, to adjust the CNTL simulation. The principle is that the CDF of the CNTL simulation should be the same as that of the observed data to adjust the total amount of precipitation. Although this bias correction has usually been used for long-term daily precipitation, we assumed that hourly data within a short-term period (approximately two weeks) could be applied as limited data were available. Note that in general, it is difficult to use the statistical bias correction to adjust the timing at which peak precipitation intensities are observed. Furthermore, the bias correction of the PGW simulation was conducted in two steps. We first conducted the statistical bias correction to which the parameters in the CNLT bias correction were simply applied and then further modified the corrected data using the linear relationship between the CNTL and PGW simulations to consider the trend of future climate change based on a similar idea by Pierce et al. [45] (Figure A2c). Figure A3 shows the accumulated amounts of precipitation among the observed and two bias-corrected precipitations for the CNTL simulation during the two weeks. The Linear BC distribution of the accumulated precipitation was approximately equivalent in amount to the observed distribution. However, the Statistical BC distribution was 10% higher than the observation. This suggests that the observed distribution could be underestimated because the kriging interpolation was performed with sparser data points than the grid spacing of the CNTL simulation, especially in the east portion of the TR watershed. Figure A2. Schematic concepts of bias correction (BC). (a) Linear BC between the CNTL simulation and observation throughs steps 1 to 3 on data arrangement in ascending order; (b) Statistical BC through the same CDF searching between CTRL simulation and observation; and (c) combination of Linear BC and Statistical BC to obtain the future trend for PGW simulation. Note that Vmod = model data, Vmdy = modified data, f = future, and α, β = coefficients.