Title Large-eddy simulation of turbulent winds during the Fukushima Daiichi Nuclear Power Plant accident by coupling with a meso-scale meteorological simulation model

A significant amount of radioactive material was accidentally discharged into the atmosphere from the Fukushima Dai-ichi Nuclear Power Plant from 12 March 2011, which produced high contaminated areas over a wide region in Japan. In conducting regional-scale atmospheric dispersion simulations, the computerbased nuclear emergency response system WSPEEDI-II developed by Japan Atomic Energy Agency was used. Because this system is driven by a meso-scale meteorological (MM) model, it is difficult to reproduce smallscale wind fluctuations due to the effects of local terrain variability and buildings within a nuclear facility that are not explicitly represented in MM models. In this study, we propose a computational approach to couple an LES-based CFD model with a MM model for detailed simulations of turbulent winds with buoyancy effects under real meteorological conditions using turbulent inflow technique. Compared to the simple measurement data, especially, the 10 min averaged wind directions of the LES differ by more than 30 degrees during some period of time. However, distribution patterns of wind speeds, directions, and potential temperature are similar to the MM data. This implies that our coupling technique has potential performance to provide detailed data on contaminated area in the nuclear accidents.


Introduction
A significant amount of radioactive material was accidentally discharged into the atmosphere from the Fukushima Dai-ichi Nuclear Power Plant (FDNPP) from 12 March 2011.In this nuclear accident, a computer-based nuclear emergency response system, Worldwide version of System for Prediction of Environmental Emergency Dose Information (WSPEEDI-II) developed by Japan Atomic Energy Agency was used to conduct atmospheric dispersion simulations from regional to hemispheric scales (Katata et al., 2012).This system consists of a meso-scale meteorological (MM) model and Lagrangian particle dispersion model, which can provide near real-time predictions of mean air concentrations and mean surface deposition of radionuclides.However, it is difficult to estimate contaminated areas in a local-scale where turbulent motions are dominant by the influence of local terrain variability and roughness elements such as individual buildings and trees within a nuclear facility that are not explicitly resolved in MM simulation models.
For simulating wind flows and plume dispersion in a localscale, a computational fluid dynamics (CFD) technique is commonly used.In CFD models, buildings and structures can be explicitly represented at high resolutions.In particular, the CFD simulations using large-eddy simulation (LES) are effective to capture complex behaviors of impinging, separating, and circulating flows around a bluff body.Therefore, an approach to couple an LES-based CFD model with a MM model is expected to have a potential of becoming an effective tool to provide detailed information on turbulent flows and plume dispersion in a local-scale under real meteorological conditions.For example, in the approach by Wys-Published by Copernicus Publications. zogrodzki et al. (2012), the MM outputs such as pressure, wind velocity, and potential temperature calculated based on RANS simulations were directly imposed at the inflow boundaries of the LES-based CFD model.Although they provided reasonable results in comparison to field experimental data, the inflow data did not include high-frequency turbulent fluctuations appropriate to drive LES-based CFD models.Michioka et al. (2013) coupled a micro-scale LESbased CFD model with a meso-scale LES-based model and investigated spatial distributions of plume concentrations in a residential area.However, they assumed a neutral atmospheric stability condition in the micro-scale model.The application of their approach is limited to a case buildinginduced turbulence is dominant.
In Japan, most nuclear facilities are located at coastal complex terrains.In this case, the assumption of neutral stability is less valid and thermal stability effects should be considered.Therefore, it is important to generate thermallystratified boundary layers with small-scale wind fluctuations in order to more faithfully reproduce meteorological conditions in the LES model from the MM outputs.In this study, we propose a calculation approach to simulation of thermally-stratified boundary layer flows by coupling an LES-based CFD model with a MM model and examine the effectiveness of the approach in comparison to measurement and MM simulation data.

CFD model
The model used for local-scale detailed simulations is the LOHDIM-LES developed by Nakayama et al. (2014).The basic equations are the filtered continuity, Navier-Stokes, and temperature transport equations under the Boussinesq approximation.The subgrid scale (SGS) turbulent effect is represented by the standard Smagorinsky model (1963) with a constant value of 0.1.The SGS scalar flux is also parameterized by an eddy viscosity model and the turbulent Prandtl number is set to a constant value 0.71.The turbulent effects of local terrain and plant canopy are represented by the external force term and are incorporated into the Navier-Stokes equation.The terrain effects are represented by immersed boundary method proposed by Goldstein et al. (1993) as follows; where m and n are negative constants.The stability limit is given by t < where k is a constant of order 1.The plant canopy effects are expressed as follows; where C d is a drag coefficient with a constant value of 0.2.U is a wind speed.a (z) is a plant area density and is determined by the forest leaf area index, thus, LAI = h 0 a (z) dz, where h is the canopy height.
Nudging terms for wind velocity and temperature fields are incorporated into the Navier-Stokes equation in order to maintain the mean structure of the MM model in the LES computational domain and can be expressed as the follows, respectively; (3) where u MM, i and θ MM are the wind velocity and potential temperature of the MM model, respectively.C nud is a spatially dependent nudging constant (see details in Sect.3.2).
The coupling algorithm of the velocity and pressure fields is based on the marker-and-cell method with the secondorder Adams-Bashforth scheme for time integration.The Poisson equation is solved by the successive over-relaxation method.For the spatial discretization in the basic equations, a second-order accurate central difference scheme is used.

Generation of turbulent inflows from MM outputs
Mean wind directions are not always constant and often vary due to a change of weather conditions.Therefore, first, we proposed the treatment of the inflow boundary conditions in order to automatically input wind velocity data obtained by the MM model into the LES model under neutral stability conditions, depending on mean wind directions in the meteorological field (Nakayama et al., 2012(Nakayama et al., , 2015)).Figure 1 shows a schematic diagram of the treatment of inflow boundary conditions in the LES model.First, mean wind directions in the meteorological field are estimated and vertical planes of inflow boundaries are determined automatically depending on them.For example, when mean wind directions α of a MM model range from −315 to 0 • or from 0 to 45 • , vertical boundary planes in the west, north, and east sides are automatically set to inflow boundaries and that in the south side is set to an outflow boundary.For the inflow boundaries, the MM outputs linearly interpolated on the spatial resolution of the LES domain.At the same time, only for the north vertical boundary plane, turbulent fluctuations generated by a turbulent inflow technique are added to the mean inflow as shown Fig. 1a.
In order to drive an LES model, a recycling method proposed by Kataoka and Mizuno (2002) is adopted.In this turbulent inflow technique, only fluctuating components are extracted at a recycling station and recycled back to the inlet boundary.The formulation is as follows: Where u, v, and w are the wind components of the streamwise (x), spanwise (y), and vertical (z) directions, respectively.The suffixes of inlt and recy indicate the instantaneous wind velocity at the inlet and the instantaneous wind velocity at the recycle station, respectively.The recycle station is set at 1 km downstream position from the main inflow boundary.
[u], [v], and [w] are horizontally averaged winds over the driver domain and ϕ (z) is a damping function.u mean , v mean , and w mean are given by mean wind velocities for each component obtained by a MM model.In cases of α ranging from 45 to 135 • , from 135 to 225 • , and from 225 to 315 • , the method to determine inflow and outflow boundaries depending on α is shown in Fig. 1b-d.
In order to produce thermally-stratified boundary layers, vertical profiles of potential temperature data obtained by a MM model are imposed at a distance of 1 km inward from each horizontal boundary.As well as the case for a wind velocity field, vertical planes of inlet and outlet boundaries are automatically determined depending on mean wind directions in the meteorological field.

Meso-scale meteorological simulation
The meso-scale meteorological simulation model used here is the Weather Research and Forecasting (WRF) model, the Advanced Research WRF Version 3.3.1 (Skamarock et al., 2008) to provide the input data for the LES model.We use a nesting capability to resolve the FDNPP region at a fine grid spacing by setting two-way nested, four computational domains (with the top being at the 50 hPa level).The four domains cover areas of 2025 km by 2025 km at 4.5 km grid, 720 km by 720 km at 1.5 km grid, 150 km by 180 km at 500 m grid, and 50 km by 50 km at 100 m grid, respectively (Fig. 2a-d).The number of vertical levels is 53, with 15 levels in the lowest 1 km depth.
The terrain data used are the global 30 s data (GTOPO30) from the US Geological Survey for the outer 2 domains and the 50 m mesh digital elevation model (DEM) dataset by the Geographical Survey Institute (GSI) of Japan for the inner 2 domains.The land-use/land-cover information is obtained from the 100 mesh dataset from the Ministry of Land, Infrastructure, Transport and Tourism of Japan.
To determine the initial and boundary conditions, we use 6-hourly Mesoscale Analysis (MANAL) data of Japan Meteorological Agency (JMA), 6-hourly Final Analysis data of the US National Centers for Environmental Prediction (NCEP FNL), and daily Merged Sea Surface Temperature (MGDSST) analyses of JMA.The times of the 6hourly MANAL and NCEP FNL are 00:00, 06:00, 12:00, and 18:00 UTC.The horizontal resolutions of MANAL and MGDSST are 10 km and 0.25 degree, respectively.Full physics processes are included in the present simulation www.adv-sci-res.net/12/127/2015/ in order to reproduce real meteorological phenomena.A physics parameterization closely relevant to the simulation of wind fields is a PBL mixing parameterization.We choose a Mellor-Yamada Level 2.5 scheme of Janjic (2002) in which mixing is done vertically between the adjacent vertical levels.A single-moment, 6-category water-and ice-phase microphysics of Hong and Lim (2006) is employed for cloud and precipitation processes in all the domains.
The case studied is the FDNPP accident on 12 March 2011.In order to simulate wind fields for this event, the time period of the WRF simulation is from 00:00 UTC 11 March 2011 to 00:00 UTC 13 March 2011.The WRF outputs during 05:00 UTC 12 March and 06:00 UTC 12 March are used for the LES model.The simulated outputs of the innermost domain at 1 min interval are used as the inputs of the LES model.

LES computational conditions
The size of the computational domain is 11.0 km by 11.0 km in the horizontal directions with the depth of 1.6 km (Fig. 2e).The total mesh number is 550 by 550 by 94 nodes.The grid spacing is 20 m in the horizontal directions and 2.5-64 m stretched in the vertical direction based on an orthogonal grid system.In the previous study, we conducted LESs of urban boundary layer flows in the urban central district by coupling with a MM model and showed reasonable results in comparison to the field experimental data of vertical profiles of wind speeds and directions (Nakayama et al., 2015).Those calculation conditions such grid resolution and computational domain were almost the same as the present ones.Therefore, it is considered that the present model set-up is reasonable to reproduce basic characteristics of the meteorological conditions in the LES model.For a wind velocity field, the inlet and spanwise boundaries are determined by the WRF wind velocity data (with 1 min interval and 100 m resolution) linearly interpolated on the spatial resolution of the LES domain with 1 min interval.At the outlet boundary, a free-slip condition is applied for each component of wind velocity.At the upper boundary, a free-slip condition for the horizontal velocity components and zero-speed condition for the vertical velocity component is imposed.For a potential temperature field, the bottom, ground surfaces, and spanwise boundaries are determined by the WRF potential temperature data (with 1 min interval and 100 m resolution) linearly interpolated on the spatial resolution of the LES domain with 1 min interval.At the outlet, a free-slip condition is imposed.
Focusing on the ground surface of the study site shown in Fig. 3a, it is found that the FDNPP is located along the coast and many forest canopies are densely situated over the land.Ground surface geometries are represented using the 50 m mesh DEM of GSI linearly interpolated on the spatial resolution of the LES model.Buffer zones with a length of 1.0 km is set in only land area and roughness blocks are placed in order to represent roughened ground surface as shown in   The nudging constant is set to rapidly vary from 0.0 to 0.01 s −1 across 750 m height using a hyperbolic tangent function for only flow field and decrease from the lateral boundaries toward the inner part using a ten-grid-point buffer zone for both flow and temperature fields.In real meteorological fields, mean wind directions are often largely different between upper and lower parts of boundary layers, which often induces numerical instabilities.Therefore, in case spatially-averaged wind directions at heights greater than 750 m height differ from those at heights less than the height by 30 degrees at the main inlet boundary, the values for each component of wind velocities are set to those at 750 m height.The time step interval is 0.05 s.The simulation period is from 05:00 UTC 12 March 2011 to 06:00 UTC 12 March 2011.2011), all monitoring posts at the FDNPP did not work due to blackout caused by the severe earthquake.Therefore, wind speeds, wind directions, and radiation dose were measured by monitoring cars.Although these simple measurement data are not appropriate to evaluate the model performance, we use the data for a comparison.The measurement data of wind speeds vary within the range from 2.7 to 3.5 m s −1 .The WRF data considerably exceed the measurement data due to the difference of the measurement height.The LES 10 min values are found to be comparable to the measurement data although the instantaneous values highly fluctuate.The measurement data of wind directions vary within the range from South-southeast to South directions.The WRF wind directions vary around Westsouthwest direction.The LES instantaneous values highly fluctuate as well as those of wind speeds.The LES 10 min values differ by more than 30 degrees from the measurement data during some period of time.
Although the difference between the LES and measurement data are observed at a ground-level during some period of time, it is successful in generally maintaining the simulated meteorological fields as the basic flow with buoyancy effects in the LES domain.

Conclusions
We proposed a calculation approach to couple the LES-based CFD model with the MM model for detailed simulations of turbulent winds with buoyancy effects under real meteorological conditions using turbulent inflow technique and examined the effectiveness in comparison to the simple measurement and MM simulation data.Inflow boundary conditions were set to automatically input wind velocity data obtained by the MM model into the LES model, depending on mean wind directions in the meteorological field.In generating thermally-stratified boundary layers from the MM outputs, first, small-scale wind fluctuations were generated at the main inlet boundary by a recycling technique.Then, the potential temperature profiles obtained by the MM model were imposed at the recycle station.
Compared to vertical profiles of the MM simulation data, it is seen that the LES wind speeds, directions, and potential temperature generally fluctuate around the MM data al-though locally rapid variations of wind speeds and directions are not reproduced well.The 10 min averaged wind directions of the LES differ by more than 30 degrees from those of the measurement data during some period of time although the averaged wind speeds of the LES are in good agreements with them.Important issues still remain in accurately simulating local-scale turbulent winds at a ground-level under real meteorological conditions.However, our approach is successful in generally reproducing thermally-stratified boundary layer flows with turbulent fluctuations in the LES model.It can be concluded that our coupling technique has potential performance to provide detailed data on contaminated area in the nuclear accidents.

Figure 1 .
Figure 1.Schematic diagram of the treatment of inflow boundary conditions in the CFD model depending on mean wind directions in the meteorological field.The solid and dotted lines are computational regions to drive turbulent winds and simulate thermallystratified boundary layers, respectively.The open and filled arrows indicate input of wind velocities and potential temperature obtained by a MM model, respectively.

Figure 2 .
Figure 2. Computational areas of the WRF and LES models.The WRF is configured with for nested domains covering areas of (a) 2025 km × 2025 km at 4.5 km grid, (b) 720 km × 720 km at 1.5 km grid, (c) 150 km × 180 km at 500 m grid, and (d) 50 km × 50 km at 100 m grid.The LES model covers an area of (e) 11 km × 11 km at 20m grid.

Figure 3 .
Figure 3. (a) The photograph reproduced by Google T M earth graphic.(b) Configuration of the FDNPP in the LES model.The buffer zone with 1.0 km is set up from each boundary.In the land area of this buffer zone, roughness blocks are placed.Green area indicates forest canopies.

Figure 4 .
Figure 4. Instantaneous fields of (a) wind speed and (b) potential temperature.

Fig. 3b .
Fig. 3b.The forest canopy is arranged based on the mixed forest in USGS 24-category land use.The LAI and canopy height are set to 4.0 and 12.0 m, respectively.The nudging constant is set to rapidly vary from 0.0 to 0.01 s −1 across 750 m height using a hyperbolic tangent function for only flow field and decrease from the lateral boundaries toward the inner part using a ten-grid-point buffer zone for both flow and temperature fields.In real meteorological fields, mean wind directions are often largely different between upper and lower parts of boundary layers, which often induces numerical instabilities.Therefore, in case spatially-averaged wind directions at heights greater than 750 m height differ from those at heights less than the height by 30 degrees at the main inlet boundary, the values for each component of wind velocities are set to those at 750 m height.The time step interval is 0.05 s.The simulation period is from 05:00 UTC 12 March 2011 to 06:00 UTC 12 March 2011.

Figure 5 .
Figure 5. Vertical profiles of wind speed, wind direction, and potential temperature obtained at (a) Main gate, (b) MP4, and (c) MP8 at 06:00 UTC 12 March 2011.The locations of Main gate, MP4, and MP8 are shown in Fig. 3.
Figure 4 shows instantaneous fields of (a) wind speed and (b) potential temperature.It is found that small-scale fluctuations in both wind and temperature fields are reproduced by the turbulent inflow technique.Figure 5 compares the LES results with vertical profiles of wind speeds, wind directions, and potential temperature of the WRF model obtained at (a) Main gate, (b) MP4, and (c) MP8 at 06:00 UTC 12 March 2011.The LES wind speeds and directions are found to be generally distributed along the WRF data and considerably fluctuate up to 100 m height at each point due to the turbulent effects by local terrain variability and forest canopy.However, at MP4, locally rapid variations of wind speeds and directions around 750 m are not captured in the LES model.The vertical profiles of potential temperature are similar to those of the WRF model.Figure6compares the LES results with the simple measurement and WRF data of time series of wind speeds and directions.The measurement and LES data are obtained at the ground-level and the WRF data are obtained at a height of

Figure 6 .
Figure 6.Time series of (a) wind speeds and (b) wind directions obtained at Main gate from 05:00 to 06:00 UTC 12 March 2011.The observed and LES data are obtained at the ground-level.The WRF data are obtained at the height of 50 m.The location of Main gate is shown in Fig. 3.