Relative Risk Assessment for Hypothetical Radioactivity Emission at a Snow Climate Site

A set of atmospheric dispersion-deposition model (ADM) integrations were conducted with a hypothetical emission of radioactive materials consisting of 137 Cs, 131 I, and 134 Cs from the Tomari Nuclear Power Plant (NPP) in Hokkaido, Japan, which is a snow climate site. Each integration was driven by Japan Meteorological Agency’s (JMA’s) mesoscale model (MSM) analysis data with 5-km horizontal resolution. The initial conditions were those on each day from January 2010 to December 2016, and the integration period was at most four days. The target was the area within 30 km of the plant. Extending a unit-mass emission concept, the measure of the relative risk is the probability of exceeding the threshold of the maximum effective dose rate based only on exposure from groundshine. Considering that the measure increased monotonically with the ratio of the total emission amount to the threshold, we evaluated the probabilistic risk with its median. The results suggested that the risk was higher in the eastern part of the target area owing to the prevailing westerly flow. The frequent snowfall in winter drags radioactive materials down in the target region, even under an active turbulent condition with strong vertical shear. The composite analysis of wind direction averaged over the target area revealed that the risk was high in the leeside but that mountains effectively blocked the inflow of the radioactive materials. The results were insen sitive to wet deposition parameterization. The risk was reduced when we replaced the emission altitude with a higher one than the standard setting. The snow shielding effect was negligible on the short-term radioactivity just after the emission but was substantial on the seasonal change in radioactivity.


Introduction
The massive release of radioactivity in the atmosphere caused by the accident at the Fukushima Daiichi Nuclear Power Plant (NPP) in March 2011 has been evaluated in many publications. Owing to the difficulty of bottom-up estimates of the emission from the nuclear reaction inside the plant, an atmospheric dispersion-deposition model (ADM) combined with radioactivity monitoring data was used to assess the time-evolution of various radioactive materials emitted from the plant inversely . Since the preliminary estimation by Chino et al. (2011), modeling studies have steadily updated the emission information with an accurate simulation of the surface deposition of radioactivity (Terada et al. 2012;Saunier et al. 2013;Winiarek et al. 2014;Katata et al. 2015; and a comprehensive review by Mathieu et al. 2018). Compared with earlier studies (Yasunari et al. 2011;Stohl et al. 2012), recent studies have input finerresolution meteorological data into more sophisticated schemes, particularly for wet deposition (Morino et al. 2013;Katata et al. 2015). Wet deposition is a crucial process in surface deposition, but the corresponding physics and chemical reactions are so complex that the wet deposition estimated with parameterization is still uncertain (Science Council of Japan 2014; Quérel et al. 2015). Moreover, because of complex topography, the observed surface deposition was entirely different from the estimation with a Gaussian plume model and flat topography (Leelössy et al. 2011). The surface deposition was influenced by local wind and precipitation patterns (Sekiyama et al. 2015, and following studies ongoing), deposition processes along the wind side of small-scale mountains, and fog or cloud deposition, as indicated by Hososhima and Kaneyasu (2015).
Based on the evacuation response to the accident of the Fukushima Daiichi NPP, the Japanese Nuclear Regulation Authority (NRA) was established and enforced nuclear emergency response guidelines, which follow the safety standards determined by the International Atomic Energy Agency (IAEA). According to the guidelines as of October 2018, at the emergency action level "general emergency", residents in the Precautionary Action Zone (PAZ), within 5 km of the NPP, should start evacuation. Residents in the Urgent Protective action planning Zone (UPZ), between 5 and 30 km from the plant, should implement protective actions such as seeking nearby shelter. After the release of radioactive materials, the Operational Intervention Level (OIL) is set on the basis of environmental radioactivity monitoring that measures the in situ effective dose rate from groundshine. The guidelines also recommend evacuation within hours in the OIL-1 area exceeding 500 μSv h −1 and a temporary relocation with an intake restriction in the OIL-2 area exceeding 20 μSv h −1 . Radioactive nuclide concentrations in food and drinks are detected in the area exceeding 0.5 μSv h −1 . Each local government where NPPs are located publishes detailed nuclear emergency planning, including annual evacuation drills, deployment of environmental monitoring posts, prescription of stable iodine, and designated evacuation routes.
The probabilistic risk assessment (PRA) of potential accidents based on an analysis considering radio active material transport (Rentai 2011) can be performed with an ADM. Studies related to the PRA include a hypothetical accident during the recovery of a Russian submarine (Brown et al. 2016), public member dose assessment of an Iranian NPP (Zali et al. 2017), and meteorological uncertainty of ADM results (Sørensen et al. 2014). A simple PRA for 16 NPPs in Japan was published by the NRA (Nuclear Regulation Agency 2012), but the PRA is not always reasonable for detailed planning including UPZ because meteorological data are sampled from a one-year observation at an NPP site, at which the surface wind is generally different from the flow that determines the advection direction of aerosols. Moreover, the PRA did not implement the wind block effect by orography or some snow effects. A common problem with these PRAs is the arbitrariness of the emission amount and radioactive nuclide types in the emission. For example, an inventory from the Fukushima Daiichi NPP was used as a severe accident in the NRA estimation, but this was simply one case. According to the emergency plan, an evacuation is decided with the OIL criteria, determined by exposure from groundshine. In general, these criteria are tunable, such as when the acceptable level of radiation exposure is changed. That is, a large emission with a severe criterion would result in a high-probability risk, and vice versa. An absolute risk evaluation is then impossible under these arbitrary parameters, and it is better to apply a unit-mass emission to the PRA.
In this study, we extend a unit-mass emission concept for PRA to present a relative evaluation of risk, and we apply this method to the Tomari NPP (Hokkaido Electric Power Co., Ltd.), located on the west coast of the Shakotan Peninsula, Hokkaido, Japan (Fig. 1). This NPP is out of operation, and a safety review is currently being conducted for a restart as of October 2018. A set of ADM integrations driven by meteorological data with 5-km resolution are performed to obtain detailed geographical distributions of deposition and atmospheric burden; the meteorological data resolution expectedly brings about not a fine result inside of PAZ but an overall result in the UPZ. The effective dose rate based on exposure from groundshine is computed from radioactivity on the ground surface with a simple formula described in Section 2. In this paper, we only examine the exposure from groundshine, on which the government determines the short-term nuclear emergency response, although a human dose assessment should include ingestion from contaminated foodstuffs and inhalation of radionuclides. Internal exposure could be formally estimated, but we would require numerous hypotheses on what residents or tourists would ingest and how they would evacuate, because the human dose strongly depends on the human activity and individual predisposition, including age and gender. Hence, this paper excludes the estimation of internal exposure, which is covered by radiology. Moreover, the snow effects are implemented in the effective dose rate estimation from groundshine because the snow season lasts for more than five months around this NPP. Section 2 also presents the meteorological data, the ADM with a treatment of snow effects, and the details of the experiments. Section 3 introduces our PRA method to provide a relative evaluation of risk with the mathematical statistics in Appendix. Section 4 briefly documents the characteristics of the local climate. Section 5 presents the geographical distribution of the relative risk and its dependence on the season, precipitation, and wind. The effects of atmospheric turbulence, snow, and emission height are also described. A risk assessment with a different threshold is provided in the Supplements. Section 6 provides a discussion related to evacuation planning and using the ADM beyond the target area. Section 7 concludes this paper.

Data and model experiments
The ADM that we used in this paper is a Lagrangian particle model including horizontal and vertical diffusion, advection by flow, gravitational settling, dry and wet deposition, and radioactive decay. This model was validated, and most of the optimal parameters were determined with a hindcast simulation for the Fukushima Daiichi NPP accident in March 2011 [see Kajino et al. (2016) for more details on the schemes and validation]. This study uses the same model framework as that of Kajino et al. (2016), although they only predicted 137 Cs and used a common scavenging parameter for both rain and snow. We included 131 I and 134 Cs in addition to 137 Cs as Lagrangian particles discharged into the atmosphere. The half-life periods are shown in Table 1. Moreover, we distinguished snow and rain on the basis of relative humidity, RH (%), and air temperature, T (°C), just above the surface. Snow falls if RH £ 92.5 -7.5T, and rain falls otherwise (Matsuo et al. 1981). We replaced the original scheme of wet deposition by Kajino et al. (2016) with that by Scott (1982) in the standard setting of this paper. The scavenging rate Λ (h −1 ) is calculated as a function of precipitation rate (P, mm h −1 ) as follows: Sensitivity experiments were performed with another parameterization with scavenging coefficient formulae: Λ = 0.30 P 0.79 for rainfall (Baklanov and Sørensen 2001), (2) Λ = 0.115 for snowfall (Kyrö et al. 2009), both assuming that the radioactive nuclei have a uniform size of 0.3 µm.
We performed ADM experiments with the hypothetical emission of radioactivity from the Tomari NPP located at 140.513°E and 43.036°N. The ADM in a single integration runs for four days or until all the particles have fallen or have been left from the calculation domain, and the radioactive materials are constantly released 20 m above the NPP during the first 24 h in the standard setting. The release rate of Lagrangian particles is 120 h −1 . The emission ratio of 137 Cs, aerosol-form 131 I, gaseous-form 131 I, and 134 Cs is 1:5:5:1, similar to an inventory for the Fukushima Daiichi NPP accident (Katata et al. 2015). We excluded gaseous-form 133 Xe, of which the total amount released was huge in the Fukushima Daiichi NPP accident, because it hardly fell to the ground and did not affect the monitoring of exposure from groundshine substantially. We also excluded other minor nuclides, such as 90 Sr and 132 Te, regardless of their lifetime. A topography with a resolution of 5 km was used in the ADM (Fig. 1a). As shown in Fig. 1a, the topographic data include Mts. Yobetsu and Inakura-Ishiyama to the north, Mts. Raiden and Niseko to the south, and Mt. Sankaku to the east of the Tomari NPP.
A set of ADM experiments were performed with initial dates every day from 1 January 2010 to 31 December 2016. Because the e-folding times of the autocorrelation function for precipitation, zonal wind, and meridional wind are 6, 27, and 18 h, respectively (not shown), the result of each integration can be regarded as independent. The ADM is driven with the grid point value mesoscale model (MSM) analysis data of the Japan Meteorological Agency (JMA), which has been calculated with the JMA's operational system developed according to Honda (2005). The horizontal grid resolution is approximately 5 km and the time interval is 3 h. The horizontal wind vector, air temperature, humidity, and precipitation variables in the domain including Hokkaido Island, as shown in Fig. 1b, are imposed on the ADM model. The output time interval is 1 h. The gridmesh for the output of the accumulated total deposition (units: Bq m −2 ) of 137 Cs, 131 I, and 134 Cs is the plane-polar grid cells with the origin at the NPP (Fig. 1a). We restricted the results to the PAZ and UPZ to maintain statistical stability. The radioactive materials were retained on the surface above the ocean to the west of the NPP (Fig. 1a), only because we wanted to see the advection effect of the easterly. To convert surface deposition data, D (Bq m −2 ), to the effective dose rate based on exposure only from groundshine, S (μSv h −1 ), we used simple formulae (Table 1 We implemented the snow shielding effect; When a radioactive snow particle is accumulated at a depth, radioactivity is reduced by 1 % by passing through 1 cm of the snowpack above (Kajino et al. 2016). The accumulation, packing, and melting of snow are numerically calculated with a simple snow depth model given by Inatsu et al. (2016). This model parameterizes snow accumulation by assuming that the fresh snow density is 0.1 g cm −3 , snow packs at a rate of 2 % day −1 , and snow melts above the freezing point at a rate of 1 cm (°C day) −1 for snow depths greater than 20 cm and at 2 cm (°C day) −1 otherwise. The removal of radioactivity by snow drift, interstitial waters penetrating the snowpack, and water immersion at the snowpack base was not considered. We also roughly implemented the snow exclusion effect: Snow particles selectively incorporate aerosol forms but exclude gaseous forms in the wet deposition process. The standard setting completely switches off the effect of gaseous-form 131 I to the groundshine estimation in snowy days, by neglecting the dry deposition effect.

Relative risk assessment
Here, we extend a unit-mass emission concept for PRA to evaluate the relative risk. First, for a given emission amount, E 0 , a single ADM integration and simple conversion formulae (Eq. 3) provide a time series of effective dose rates for a particular grid point, corresponding to the radioactivity measured by a monitoring post in the evacuation planning. For example, if the maximum value at the point in an integration exceeds 20 μSv h −1 , the point satisfies the OIL-2 criterion. Following this concept, the number of integrations, n, where the maximum effective dose rate exceeds a criterion C 0 , is counted at every grid point shown in Fig. 1a. The probability of risk is then evaluated as where the total number of integrations is N = 2,557. The conditional probability in the composite analysis, by sampling integrations to meet a condition, can be defined similarly. We considered a set of maximum effective dose rates at a fixed grid point in each integration with a descending order, The risk probability is The emission amount is α-fold, and then the same risk probability is obtained for the α-fold criterion, because αC n 0 > αC 0 ³ αC n 0 +1 . Extending this discussion, the α-fold emission amount is equivalent to the α −1 -fold criterion. Hence, we naturally introduced a single parameter of the ratio of emission amount to the criterion, hereafter called the E/C ratio, and only the emission amount of 137 Cs was used because the ratio among 137 Cs, 131 I, and 134 Cs is fixed at 1:10:1.
The probability can be calculated for a given emission amount and a given criterion of effective dose rate and then represented as a function of the E/C ratio and grid point. We call the line plot for a grid point the risk curve because the concept in the original use of this terminology in insurance is similar. Figure 2 displays the risk curves for each grid point in the target area based on the ADM integrations. A higher probability with a fixed E/C ratio and a lower E/C ratio with a fixed probability both indicate a higher risk; for example, the upper, leftward curve in Fig. 2 shows a higher risk. The risk curve for a grid point in the PAZ shows the highest risk. The risk in the eastern quadrant of the UPZ is much higher than that in the western quadrant. Because the risk curve is not always represented as the same function form, the comparison between two risk curves is generally difficult unless the probability is fixed. However, the result in Fig. 2 suggests that the Z-value associated with the probability for each grid point has a similar shape as a function of the E/C ratio: The curve is steep on the negative side of Z, whereas it is gentle on the positive side. This notion is supported by a mathematical consideration of random sampling from the analytic solution of a two-dimensional linear advection-diffusion equation (Appendix). The aim of the method of relative risk evaluation in this paper is to show the E/C ratio at a 50 % probability, that is, the median, on a map. To evaluate the relative risk at a higher or a lower probability, the fitting function in Appendix might be helpful, but we provided the risk evaluation with the E/C ratio at a 5 % probability in the Supplements for the readers' convenience.

Local climate
Here, we briefly describe the climate around the target area (Fig. 3). The low-level westerly prevails in this area throughout the year related to the subtropical jet stream above the planetary boundary layer (PBL).
The summer climatological wind is southwesterly in the northwestern edge of the Bonin high. The annual precipitation is around 1000 mm near the NPP, but precipitation is more than 2000 mm at Mts. Yobetsu and Inakura-Ishiyama in the Shakotan Peninsula and at Mts. Raiden and Niseko (Fig. 1a). Kutchan, the capital of the Shiribeshi Subprefecture (Fig. 1a), exhibits heavy snowfall with an annual maximum snow depth near 200 cm, mainly due to cloud bands brought in by the westerly in the cold season. In contrast, snowfall is lighter along the coastal line, with an annual maximum snow depth of only 73 cm at Suttsu (Fig. 1a).

Overall features and seasonal variation
The results of the ADM integrations are presented as relative risk for radioactivity with the median of the E/C ratio as the index. The annual mean risk is distributed anisotropically from the NPP where radioactive materials were hypothetically released (Fig. 4a). The high-risk area spreads to the east because the radioactive materials emitted from the NPP are advected by westerly wind (Fig. 3). Despite many possible weather patterns that produce a southerly, the risk is low in the leeside of the Shakotan Peninsula. Similarly, even though a northerly pattern would mainly send (shading) annual snowfall estimated with an assumption of rain snow discrimination following Matsuo et al. (1981), and (vectors) annual mean wind vector at 950 hPa, all based on MSM analysis data that we provided in the ADM integrations as the forcing. The contour interval is 100 mm with the 1,500, 2,000, and 2,500 mm lines thickened; the shading is as per the reference in the right; and the reference vector is 5 m s −1 in the bottom right corner. particles south, the high-risk area is limited around the NPP. These results are explained by Mts. Yobetsu and Inakura-Ishiyama in the north and Mts. Raiden and Niseko in the south (Fig. 1a), blocking the inflow of the radioactive materials. The no-deposition days show a similar horizontal distribution with counts of less than 20 % in the eastern quadrant of the UPZ and 40 -60 % in the western quadrant (Fig. 4b). Figure 5 shows the relative risk in each season. In boreal winter (Fig. 5a), the relative risk is the highest. The distribution of the risk area is similar to the annual mean, but the median of the E/C ratio is smaller than 1 PBq (µSv h −1 ) −1 in most of the eastern quadrant. This is closely related to the fact that approximately 45 % of the emitted radioactive materials would fall in the target area (Figs. 6a, 7a). More than 80 % of the maximum effective dose rate is attributed to the wet deposition in the eastern quadrant, and dry deposition is much less effective even in the ocean area (Fig. 6a). In contrast, the relative risk is much lower in boreal summer (Fig. 5c). Only a few grid points near the NPP satisfy < 1 PBq (µSv h −1 ) −1 . This is related to less than 20 % of the emission falling in the target area with dry deposition explaining about half of the maximum effective dose rate (Figs. 6b, 7b). In contrast to other seasons, the risk area expands not only east, but also northwest. The risk is low in the north of the Shakotan Peninsula and south of the NPP. Spring and autumn show relative risks intermediate in character between those of winter and summer (Figs. 5b, d).

Effect of precipitation and wind
The overall features and seasonal variation of the relative risk distribution are attributed mainly to wet deposition by precipitation and partly to dry deposition by settlement due to atmospheric turbulence, which depends on static stability and vertical wind shear. Horizontal wind, especially in the PBL, controls the direction and speed of particles. Hence, precipitation and wind are two major factors that determine where radioactive particles move and fall.
About half of the total days in seven years experienced a daily precipitation amount exceeding 1 mm at the NPP (Table 2). For these days, the risk is higher than that in the annual mean (Fig. 4a); the median of the E/C ratio probability reduces to almost half (Fig.  8a). The high-risk area spreads over the eastern quadrant of the target area, similar to that of the annual mean. A quarter of days were categorized into heavier precipitation exceeding 5 mm day −1 . The composite for this category shows that the E/C ratio is even lower (Fig. 8b). The E/C ratio decreases much more, and the risk is higher near the NPP site for heavy precipitation of > 10 mm day −1 at the NPP (Fig. 8c). About 13 % of the total days fell into this category. These results are explained by the high sensitivity of the scavenging coefficient to precipitation intensity (Eq. 1). Although the index is now based on the precipitation at the NPP, it is highly correlated with the precipitation at many grid points in the target area (not shown). The risk difference between summer and winter is mainly caused by the precipitation frequency. The precip-     itation exceeded 1 mm day −1 on only one third of summer days ( Table 2). The atmospheric burden in summer reaches 20 % of the total emission amount in the first 24 h in the simulation (Fig. 7b) because the lower frequency of precipitation leaves more than 80 % of the radioactive material suspended in the air. In contrast, on 77 % of winter days, precipitation exceeded 1 mm (Table 2), increasing the wet deposition in the target area (Fig. 7a). The atmospheric burden and radioactivity leaving the target area in winter are much less than those in summer. The composite analysis for each wind direction is summarized in Fig. 9. The westerly blew in the PBL in one-third of the cases, pushing radioactive parti-cles toward the east (Fig. 9f). The northwesterly and southwesterly are the second and third frequent cases, respectively, but the northwesterly has a higher risk in its downstream region than the southwesterly (Figs. 9c, i). This is because the northwesterly is mainly observed in winter with heavy snowfall, whereas the southwesterly occurs in warm seasons (Fig. 9e). The southerly case is frequently observed in summer. However, the risk just north of the NPP is reduced by the mountains (Fig. 9b). The flow tends to bend along the western coast of the Shakotan Peninsula, so that the risk area spreads toward the northwest. Similarly, because Mts. Raiden and Niseko block the flow, the high-risk area extends to the southeast in the norther- ly composite (Fig. 9h). Other minor cases also show downstream deposition; most of the radioactivity would fall in the ocean in the easterly and southeasterly cases (Figs. 9a, d). Because the composite for the wind direction at 950 hPa, almost 500 m above sea level, explains the annual mean and seasonal variations of risk (Figs. 4a, 5), the advection of particles is mostly controlled by wind at this level (Fig. 3).
The number of particles that reach the upper atmosphere is affected by atmospheric turbulence. The climatological vertical diffusion coefficient is large at the bottom of the PBL in winter (Fig. 10a) because of the strong vertical wind shear in the baroclinic zone (Fig. 10b) and weak static stability in the PBL owing to the lake effect accompanied by shallow convec-tion (Fig. 10c). However, most particles do not reach the PBL and fall to the ground within several hours of their emission (Fig. 11a). In contrast, weak wind shear and stable stratification characterize the PBL in summer, resulting in a small diffusion coefficient (Fig.  10a). Despite the small diffusion coefficient, most particles stay in the PBL in summer (Fig. 11b). Hence, the effect of atmospheric turbulence is limited to the timing just after emission, and the degree of risk in the target area is controlled by wet deposition rather than atmospheric turbulence.
We finally examined how much the snow shielding and gaseous-form exclusion affect the results. Figure  12 shows the E/C ratio in December-January-February (DJF) estimated without the gaseous-form effect. The  risk estimate is almost the same as that with the effect (Fig. 5a), probably because of the low conversion rate for 131 I (Eq. 3). This comparison clearly shows the little impact of this gaseous-form exclusion effect on the risk assessment. We next thought of the snow shielding effect that radioactivity decreases by passing through the snowpack. This effect was measured during monitoring after the Fukushima Daiichi NPP accident, when unpolluted new snow covered the polluted ground (Kajino et al. 2016). Our experiment ran the ADM for only four days in each integration, so that the snow depth above the radioactive materials on the ground or in the snowpack would not be enough to reduce radioactivity. This effect would be less than 5 % of the wintertime average (Fig. 13a). However, the snow shielding would substantially decrease the radioactivity during the snow season after the emission date. For example, because the annual-maximum snow depth often exceeds 200 cm in Kutchan town center (Figs. 1a, 3; Katsuyama et al. 2017), the decay rate of the radioactivity would be noticeable in an integration where the emission happened in early winter. Figure 13b shows an example of the radioactivity decay rate relative to the initial day value at Kutchan for an emission on 1 December 2010. A simple snow depth model simulation extended until the end of May, forced with bias-corrected surface temperature and precipitation of the MSM analysis at the grid point nearest to Kutchan. The simulated snow depth is slightly overestimated, but the snow accumulation period is almost reproduced. Considering the ratio of radioactive nuclei types (Section 2) and their conversion rate ( Cs to the effective dose rate on the initial day were similar. After the half-life period of 131 I of eight days, the effective dose rate would drastically decay regardless of snow accumulation. Because the snow depth exceeds 100 cm in January in this example, the effective dose rate would decrease by more than 90 % and would remain at this low level throughout the snow season. The radioactivity would then increase again at the end of April, when the snow completely melts. The results suggest that the snow shielding effect is  negligible in the short-term risk assessment, which was the focus of this work, although it is substantial in the seasonal change in radioactivity.

Sensitivity test
A sensitivity test of the wet deposition parameterization was performed. The wet deposition parameterization in the standard experiment (Eq. 1) allows many particles to fall to the ground, compared with other size-resolved scavenging coefficient formulae (Wang et al. 2010;Zhang et al. 2013). We replaced the parameterization with empirical rain (Baklanov and Sørensen 2001) and snow (Kyrö et al. 2009) scavenging parameterization, and we founds a similar spatial distribution of relative risks with the empirical parameterization (Figs. 4a, 14a), although a high-risk region extends eastwards in the sensitivity experiment due to inefficient scavenging. Even in the heavy precipitation composite, the wet deposition parameterization is insensitive to the risk distribution in the UPZ.
A further sensitivity test of the emission height was performed. The emission height is fixed at 20 m above the ground level in the standard experiment, following the Fukushima Daiichi nuclear disaster (Katata et al. 2015). However, changing the emission height to 200 m does not alter the relative risk distribution greatly in the UPZ (Fig. 15). The risk in the PAZ and UPZ is lower in the experiment, presumably because particles at a greater height would move faster initially, and thus dry deposition is less effective above 50 m.

Discussion
This study provides some different results from the NRA simulation (NRA 2012 or link to a Japanese document at http://www.nsr.go.jp/data/000024448.pdf), which only used the surface wind at the NPP to drive a Gaussian plume model. The NRA simulation may be accurate within a few kilometers of the NPP because the particle motion close to the ground just after emission is controlled by near-surface wind for a short time, until the particles fall to the ground or are lifted several hundreds of meters. However, the wind direction in this simulation is strongly influenced by the mountains behind the NPP (Figs. 1a, 16a) and is different from the PBL flow (Figs. 2, 9e). Although the westerly is dominant in the PBL flow, the southerly and southwesterly are less frequent because of the mountains blocking the winds. The surface flow above the NPP tends to be aligned with the western coast of the Shakotan Peninsula, and thus the southeasterly is dominant instead. This observation is consistent with the risk distribution in summer (Fig. 5c), but the NRA result greatly deemphasized the advection toward the north. However, when the NRA performed the simulation for every-hour emission, their approach may have maintained a higher temporal resolution than ours with 24-h continuous emission, even though this was not discussed in the NRA's report. As shown in the surface wind rose in summer at the NPP (Fig.  16b), the prevailing wind is easterly in the morning and westerly in the evening. This is probably because the sea breeze blows from the colder land surface to the warmer sea surface in night to morning and it reversely blows in the evening. This diurnal cycle is prominent in summer, but not dominant in cold seasons due to the strong westerly flow (not shown). We did not deal with this problem directly, but it could be resolved if the daily emission were replaced with an hourly emission like the NRA's simulation with a pru-  dent treatment for the independence of integrations.
The spatial resolution and sampling number in our ADM experiment are sufficient to resolve the local features of wind, precipitation, and topography in the target area. The relative risk assessment in this study could formally extend out of the target area, but there are two obstacles that can reduce the robustness of results. First, this study uses 2880 Lagrangian particles emitted from the NPP. If the particles were to travel straight radially and isotropically, 80 particles could pass over a grid cell in the outermost grid of the target area. However, the number of particles that would pass over the boundary of a 100-km circle would be less than 20, which is not sufficient to maintain statistical stability. In other words, a locally highconcentration area out of the UPZ, called a hot spot, could appear, as in the Fukushima Daiichi NPP accident, but because the trajectory of particles is sparsely distributed far from the emission site, sufficient particles are needed to evaluate the probability. The PRA commonly suffers from this problem, which can be solved by intensive computation with enough particles and integrations. The other problem is related to the complex geography to the east of Tomari NPP, where a range of mountains is located, with Mt. Yoichi (> 1000 m) as the highest. Many small valleys and ridges may produce precipitation patterns along with the wind side of the ridges and a complex diffluence and confluence pattern. Simulation with 5-km resolution cannot resolve such small-scale phenomena. Hence, an additional high-resolution simulation, such as dynamical downscaling (Inatsu et al. 2015) or regional reanalysis (Mesinger et al. 2006), is necessary for a reliable risk assessment out of the target region for Tomari NPP, but this is beyond the scope of this paper.

Conclusion
We performed a set of ADM integrations, each of which contains four-day radioactive material transport, including dry and wet deposition with a hypothetical emission during the first 24 h of 137 Cs, 131 I, and 134 Cs from Tomari NPP, Hokkaido, Japan, which is a snow climate site. The ADM was driven by meteorological data with 5-km horizontal resolution from January 2010 to December 2016. Extending a unit-mass emission concept, the results were arranged in light of Japanese nuclear emergency planning, in which evacuation is triggered on the basis of the effective dose rate, which only considers exposure from groundshine measured by monitoring posts dispatched in an emergency. We demonstrated that the probability of the maximum effective dose rate exceeding a criterion is a monotonic function of the E/C ratio (Appendix). The relative risk is assessed with the median of the E/C ratio. The risk was high in the east of the NPP in autumn, winter, and spring, owing to the westerly flow prevailing in the PBL, whereas it expanded northwestwards as well in summer (Figs. 4a, 5). The composite and budget analysis revealed that wet deposition played a primary role in risk distribution (Figs. 7,8); the risk in winter with frequent heavy snowfall was higher than that in summer with many sunny days ( Table 2). The atmospheric turbulence contributed to lifting particles just after the emission, but the wet deposition dragged most of the particles to the ground (Figs. 7, 10, 11). The wet deposition parameterization was insensitive to the risk distribution (Fig. 14), whereas the emission at a higher level reduced the risk (Fig. 15). The snow shielding effect was negligible on the short-term radioactivity just after the emission but was substantial on the seasonal change in radioactivity (Fig. 13).

Supplements
Supplements provide the results with the E/C ratio at a 5 % probability.

Acknowledgments
We would like to thank the two anonymous reviewers who gave us helpful comments to improve this paper. We are deeply indebted to Prof. T. Kozaki and other members of the nuclear disaster prevention committee of the Hokkaido local government for discussions related to nuclear emergency planning in Hokkaido. Professors M. Mori, H. Ishikawa, S.

Appendix: Mathematical consideration of the risk curve
This appendix demonstrates that the risk curve shown in Fig. 2 has a steep slope on the low-probability side and a gentle slope on the high-probability side. First, we assume that the deposition follows a bivariate Gaussian with its center drifted by wind and is an analytic solution of the two-dimensional linear advection diffusion equation. We do not explicitly include the effect of precipitation, but the following consideration can partially reflect heavy rainfall, creating a high-concentration area following a bivariate Gaussian. The statistics in this paper are rewritten as a sampling from many sets of drifted data following bivariate Gaussianity at a fixed point. This is equivalent to sampling at a randomly moving point from a fixed Gaussian. Without loss of generality, the Gaussian can be normalized as follows: where r is the distance from the origin. The sampling point is assumed to follow another Gaussian of where (x 0 , y 0 ) = ρ(cos ϕ, sin ϕ) is the central point of sampling and σ is the standard deviation of this distribution. The probability density function (PDF) for this sampling is written as where θ is the angle of polar coordinate. Noting dr/dF = -1/rF, where I 0 is the modified Bessel function of the first kind and has the Maclaurin series of where b = log(2πX ) (< 0). Using Eq. (A5), the integral in part, and the Taylor series of exponential function, the asymptotic solution of the CDF is where ξ = aρ 2 /2. By taking the limit of b ® -¥, then Φ ® 1. Therefore, the PDF satisfies the condition that the total probability is unity. Finally, using Fisher's Z-transformation, the Z-value is explicitly written as a function of criterion X as follows: The graph of Z (X ) is displayed in Fig. A1. Because ρ is the distance from the sampling center to the origin, with a fixed sampling range σ, it represents the dominance of diffusion compared with the advection effect. Sampling at a point strongly affected by advection or near the emission point corresponds to sampling the data near the PDF peak. Figure A1 indicates that the effect of advection makes the risk curve steeper at the low-probability side.