A semi-Lagrangian advection scheme for radioactive tracers in the NCEP Regional Spectral Model ( RSM )

In this study, the non-iteration dimensional-split semi-Lagrangian (NDSL) advection scheme is applied to the National Centers for Environmental Prediction (NCEP) Regional Spectral Model (RSM) to alleviate the Gibbs phenomenon. The Gibbs phenomenon is a problem wherein negative values of positive-definite quantities (e.g., moisture and tracers) are generated by the spectral space transformation in a spectral model system. To solve this problem, the spectral prognostic specific humidity and radioactive tracer advection scheme is replaced by the NDSL advection scheme, which considers advection of tracers in a grid system without spectral space transformations. A regional version of the NDSL is developed in this study and is applied to the RSM. Idealized experiments show that the regional version of the NDSL is successful. The model runs for an actual case study suggest that the NDSL can successfully advect radioactive tracers (iodine-131 and cesium-137) without noise from the Gibbs phenomenon. The NDSL can also remove negative specific humidity values produced in spectral calculations without losing detailed features.


Introduction
The spectral method is a numerical method applied to primitive meteorological equations to achieve high-order accuracy (Robert, 1966).The spectral method has advantages over the finite difference method: the absence of truncation error in linear terms, the absence of aliasing and phase errors, the absence of pole problems in global models, and the efficient use and simple programming of semi-implicit time integration (Bourke, 1972).Applying the spectral method in a regional model can be difficult due to time-dependent lateral boundary conditions.However, several approaches for solving this problem have been successful (e.g., Tatsumi, 1986;Fulton and Schubert, 1987;Hoyer, 1987;Segami et al., 1989;Chen and Kuo, 1992).Juang and Kanamitsu (1994) presented the National Centers for Environmental Prediction (NCEP) Regional Spectral Model (RSM).The RSM has advantages in accuracy for a regional high-resolution domain.In addition, the spectral representation of the RSM is a two-dimensional perturbation method, which can eliminate the error due to reevaluation of the linear forcing from the base fields by the regional model (Juang et al., 1997).This is one of the reasons that the RSM can be easily used for long-range climate simulations.The RSM has been widely used for regional downscaling research (e.g., Kang and Hong, 2008;Kanamitsu et al., 2010;Chang and Hong, 2011;Li et al., 2012).Kang and Hong (2008) assessed impact of the land surface parameters on the regional climate circulations.Kanamitsu et al. (2010) presented a refined spectral nudging technique for regional dynamical downscaling.Chang and Hong (2011) used the RSM to produce regional future scenarios by dynamical downscaling.Li et al. (2012) showed that the fully coupled RSM and regional ocean modeling system (ROMS) can produce detailed oceanic circulations over the California coast.
Although the RSM has advantages, it also has a problem common in spectral model systems: the Gibbs phenomenon.This phenomenon is "overshooting" in the convergence of the partial sums of particular Fourier series in the neighborhood of a discontinuity of the function being expanded.In the case of spectral techniques, the Gibbs phenomenon can introduce negative values in positive-definite variables E.-C.Chang and K. Yoshimura: A semi-Lagrangian advection scheme for radioactive tracers (i.e., hydrometeors and tracers).Yoshimura (2011) presented a tracer simulation for the Fukushima Daiichi nuclear power plant accident during the massive earthquakes and tsunami on 11 March 2011 in eastern Japan by utilizing the stable isotope mode of the RSM (IsoRSM; Yoshimura et al., 2010).General features are well captured in this simulation, e.g., the tracers from Fukushima reached the metropolitan area and a significant amount of tracers that were precipitated.However, these simulations also clearly show that computational noise is produced by the Gibbs phenomenon near the emission point.This problem is quite severe in this simulation because radioactive materials are emitted from a single gridpoint source at the surface level, which creates a significant discontinuity in the wave space transformation.Thus, an advection method that does not require a spectral transformation for tracers is needed to resolve the Gibbs phenomenon.
The semi-Lagrangian method is an alternative advection scheme that can replace the spectral calculation of hydrometeors and tracers.Semi-Lagrangian advection schemes have long been preferred in numerical weather prediction because they are more accurate and efficient than traditional Eulerian schemes when large time steps are considered (Williamson, 2007).Staniforth and Côté (1991) reviewed the semi-Lagrangian literature for atmospheric models.They concluded that the semi-Lagrangian framework facilitates the incorporation of shape-preserving and monotonic schemes for moisture advection, because of the relatively small dispersion errors in the presence of discontinuities or near discontinuities.Juang (2007Juang ( , 2008) ) proposed a non-iteration dimensional-split semi-Lagrangian (NDSL) scheme, which is simple and economical compared with the conventional semi-Lagrangian method.Although the NDSL scheme has not yet been applied to the regional spectral model, some versions of the NDSL have been added to the global spectral model system, e.g., the NCEP Global Forecast System (GFS; Moorthi et al., 2001) and the Global/Regional Integrated Model system (GRIMs; Hong et al., 2013) Global Model Program (GMP).For a regional model system, the flux through the boundaries is needed to apply the mass restoration, whereas there are no boundaries for global domains.Aranami et al. (2015) applied a mass restoration scheme for limited-area models (LAMs) with semi-Lagrangian advection.As such, the boundary treatment is required to apply the NDSL advection scheme for the RSM.
The objective of this study is to remove the Gibbs phenomenon for hydrometeors and tracers in the regional spectral model by replacing the spectral tracer advection scheme with the semi-Lagrangian advection scheme in a real simulation.Detailed features of the semi-Lagrangian version of the regional spectral model and the experimental design are described in Sect. 2. Section 3 provides results from the original IsoRSM and the semi-Lagrangian version of IsoRSM for radioactive tracers and humidity fields in the Fukushima nu- clear power plant accident case study.A summary and conclusion are provided in Sect. 4.

Method
2.1 Regional spectral model for radioactive tracers Yoshimura et al. (2010) presented the IsoRSM, which includes the isotopic species for water vapor (HDO and H 18  2 O) as tracers in the latest version of the Scripps Experimental Climate Prediction Center's regional spectral model (Kanamitsu et al., 2005).In the IsoRSM, the tracers can be incorporated into raindrops or cloud particles at every integration time step.Therefore, the interactions between the tracers and precipitation processes can be considered.In contrast, a general chemical transport model uses precipitation as an external forcing.Yoshimura (2011) and Saya et al. (2013) modified IsoRSM to enable the simulation of the transport of radioactive tracers.Isotopic variables were replaced with radioactive tracers (i.e., 131 I and 137 Cs), and dry and wet deposition processes caused by gravity and precipitation processes, respectively, were introduced.Saya et al. (2013) proposed a wet deposition process wherein the deposition is proportional to the ratio of the amount of condensed water to the total amount of water, whereas a traditional method (e.g., Maryon et al., 1991) considers only the amount of condensed water.In this study, the radioactive tracer mode of IsoRSM is used as an original framework.

NDSL scheme for RSM
In this study, the NDSL advection scheme replaces the spectral prognostic calculation of the tracers in the RSM.Once each tracer field is provided from the initial field on the regular grid space, the advection of these fields is calculated by the NDSL on the grid space without any spectral transformation during the model integration.This process prevents the Gibbs phenomenon.The NDSL has two characteristics: (1) non-iteration to compute the trajectories of each tracer and (2) a dimensional-splitting method.Figure 1 shows the basic concept of two-dimensional advection for both the traditional semi-Lagrangian scheme and the NDSL scheme (from Juang, 2008).The traditional backward (forward) semi-Lagrangian scheme assumes that the arrival (departure) points are located on the regular model grid (Fig. 1a).In this case, an initial guess and iterations to compute the trajectories are required, which means finding a midpoint wind and transferring the fluid particles from the departure points to the arrival points.These iterations make the semi-Lagrangian scheme expensive and inefficient.The NDSL scheme is a central scheme, which assumes a midpoint wind at the regular model grid point at time t to find the departure point at time t − t and the arrival point at time t + t (Fig. 1b, also see Fig. 1 of Zhang and Juang, 2012).No initial guess or iteration is needed to find trajectories.Only one interpolation at the departure point and one remapping at the arrival point are needed.The wind in the central scheme can reduce numerical error better than the upstream scheme because no estimation or iteration occurs.During the advection process, the quantity of the tracer is assumed to be steady.The NDSL uses the dimensional-splitting method for 2-D advection, which simply splits the 2-D advection into a sequence of 1-D advections.Because it only needs 1-D interpolation and remapping, the method easily attains mass conservation.Another important advantage is that the 1-D method can be easily coded; thus, it is compatible with most numerical models.
For the implementation of the regional version of the NDSL in a real-data case, the boundary treatment must be considered.The global NDSL uses a cyclic boundary condition, which places a duplicated global domain at the western and eastern boundaries.This allows for the arrival or departure point to be located, even when it lies outside of the model domain.However, cyclic boundary conditions are not suitable for regional domains because the domain edge points differ in both longitude and latitude.Thus, the domain is treated as three sections: the boundary zone, the buffer zone, and inner domains.Figure 2 shows the structure of these sections in the lower-left corner of the model domain as an example; it assumes that the boundary and buffer zones are defined as three grid points each.The boundary zone (dark gray in Fig. 2a) is the area where the semi-Lagrangian advection calculation is not applied and where the values of the global base field are specified.This area is necessary to prevent the calculated departure point of tracers located outside of the model domain.The values in the inner domain are calculated entirely from the regional model.In the buffer zone (light gray in Fig. 2a), the global base field and the regional model field are combined, with the weighting determined by an inverse exponential function (Eq. 1, Fig. 2b), to smooth the gap between the boundary zone and the inner domain.
Here, W G is the weighting for the global base field, W R is the weighting for the regional model field, and k is the grid point of the buffer zone toward the inner domain.W G is defined according to Eq. (1) in the buffer zone but is specified as 1 and 0 in the boundary zone and the inner domain, respectively.Finally, the result over the model domain is defined by Eq. ( 2): where F is the final result of a particular tracer, F G is the global base field, and F R is the field calculated by the regional model.

Idealized experiment
An idealized experiment is performed to examine the feasibility of the regional version of the NDSL advection scheme.A horizontal advection scheme is examined in the 2-D domain, which has 400 east-west grid points and 400 northsouth grid points, with a 10 km resolution.A uniform wind field of 10 m s −1 in both the x and y directions is used for horizontal advection.The integration time interval ( t) is set to 10 s. Figure 3a shows the result of this horizontal advection every 1000 time steps.An ideal perturbation is imposed at the initial step, and the shape of the perturbation is maintained during the integration up to the 3000th time step.The transported disturbance by the NDSL scheme at the 3000th time step is compared to the analytic solution (Fig. 3b); these two results are nearly identical.The shaded values in Fig. 3b indicate differences between the analytic solution and the result calculated from the NDSL.The differences are less than 1 % of the maximum disturbance.The ratio of the mass change due to the NDSL to the initial mass (Eq. 3) is verified to confirm that the NDSL satisfies mass conservation.
R mass = (total mass) − (initial total mass) (initial total mass) (3) R mass ≈ 10 −15 up to the 3000th time step, which shows that the NDSL scheme satisfies mass conservation during the integration.
The idealized vertical advection experiment is performed in one dimension.The vertical dimension is the sigma coordinate, which has 100 layers from 1 (bottom) to 0 (top) with equal spacing (0.01).A uniform vertical velocity of 10 −4 sigma s −1 is prescribed.The integration time interval ( t) is 100 s. Figure 4 shows that the virtual concentration at the initial time is conserved during the transport process.The mass conservation ratio for vertical advection is R mass ≈ 10 −15 .The idealized experiments in the horizontal and vertical directions verify that the regional version of the NDSL scheme can accurately calculate the advection of these specific perturbations in the horizontal and vertical directions while conserving mass.

Experiment design
For the Fukushima case study, two experiments are performed.The first is the ORG run, which is identical to the control experiment of Saya et al. (2013), wherein tracer fields are calculated in the spectral space, as is performed in the original IsoRSM.The second is the SL run, wherein the humidity and radioactive material ( 131 I and 137 Cs) tracer fields are calculated by the NDSL scheme.All configurations except the advection method for tracers are identical in the two experiments.Figure 5 shows the experimental domain for the case study, with horizontal grid spacing of 10 km.The number of grid points is 161 (east-west) by 200 (north-south).The number of vertical layers is 28 in terrain-following sigma coordinates; the lowest and highest sigma levels are 0.995 and 0.002, respectively.The red circle in Fig. 5 indicates the emission point, which is the location of the Fukushima Daiichi nuclear power plant.The physical processes used are the relaxed Arakawa-Schubert deep convection scheme (Moorthi and Suarez, 1992), the Noah land surface model (Ek et al., 2003), the Chou radiation scheme (Chou and Suarez, 1994), and a non-local planetary boundary scheme (Hong and Pan, 1996).The model simulation is integrated from 00:00 UTC 12 March 2011 to 00:00 UTC 28 March 2011 (16 days).Atmospheric initial and lateral boundary conditions are provided by the NCEP-Department of Energy (DOE) reanalysis (Kanamitsu et al., 2002).In the case that negative values are introduced in the initial field, correction is performed in regional interpolation process by replacing negative values with zero.If negative tracer quantities are produced by the physical parameterizations, those negative values are transported to the layer above and the original values are replaced by zero.The emission rate of the radioactive tracers from Chino et al. (2011) is used in this study.To determine the size of the boundary and buffer zones (Fig. 2) for this case study, the fastest wave is assumed to be a sound wave, which moves at a speed of 300 m s −1 .The estimated longest traveling distance of this wave in one time step ( t = 40 s) is 12 km; this value is less than 2 x.Therefore, the fastest-moving tracer is confined within 2 grid points over one time step.Thus, it is impossible for tracers to enter from outside the domain when the boundary zone is larger than 2 x.In this study, the boundary zone is set to 5 x for safety.The buffer zone is the same size as the boundary zone.
When the NDSL advection is used (the SL experiment), the extra computation cost is 35 % higher with respect to the ORG run.However, the current version of the NDSL in the RSM still calculates spectral tracer advections even though the result is not used any more.It means that the computational burden with the current NDSL in the IsoRSM is purely the increased computational cost which is required for the NDSL tracer advections.In an updated release, this inefficiency will be solved.

Radioactive tracer field
Figure 6 shows column-integrated atmospheric radioactive tracer ( 137 Cs) from the ORG and the SL experiments at 12:00 UTC, 15 March 2011, when the maximum emission rates occurred.The ORG run produces a very distinct noise in the zonal and meridional directions from the emission point (Fig. 6a).Additionally, a ring-shaped signal surrounds the emission point.This signal is the pattern of highconcentration, empty values; the ring shape extends from the center of the emission to outside the domain.These signals, also known as the ringing artifact, are typical of the Gibbs phenomenon.In the Fukushima case study, the discontinuity of the tracer field is pronounced because the tracers are emitted from a single grid point, which leads to significant noise from the spectral transformation processes in the ORG experiment.Widespread distributions of tracers with small concentrations occur over the domain.However, the SL experiment does not have any noise from the Gibbs phenomenon (Fig. 6b).Compared with the ORG run, the SL experiment produces a generally similar pattern of tracers, which advect  7 shows longitudinal-vertical cross sections of tracers averaged in the meridional direction over 36.5-37.5 • N, which includes the emission point.This figure shows that the ringing noises extend above the surface in the ORG experiment (Fig. 7a).The SL experiment shows transport in the vertical and horizontal directions without computational noise.Clearly, the semi-Lagrangian advection method can calculate the transport of tracers without computational noise, whereas the spectral representation of advection produces severe errors.

Humidity field
The SL experiment applies the semi-Lagrangian method to radioactive tracers and the humidity field.The humidity field is also a positive-definite field and spectrally exhibits the Gibbs phenomenon, although the discontinuity is not as strong as it is in the radioactive tracer fields.Figure 8 shows that the specific humidity in the ORG run exhibits some negative values in the lower troposphere (Fig. 8a, b).The SL experiment simulates detailed humidity distributions that are similar to those of the ORG experiment but without any negative values (Fig. 8c, d).Negatives in the ORG run indicate that the original RSM has a systematic problem representing the positive-definite field, even though the RSM is widely used and evaluated for regional downscaling.The 16-day accumulated precipitation fields from each experiment are quite similar (Fig. 9a, b).General rainfall patterns observed in the Tropical Rainfall Measuring Mission (TRMM) Multi-satellite Precipitation Analysis (TMPA) are well captured in both experiments (Fig. 9c).The spatial correlation coefficient of precipitation between the ORG run and the TMPA is 0.616, whereas the correlation coefficient between the SL run and the TMPA is 0.622.It means that the corrected humidity field by the NDSL scheme can slightly improve precipitation or keep the simulation skill of the original IsoRSM in the rainfall simulation.When we consider that the ORG experimental set has been widely used for various downscaling research, it is possible to understand that the regional NDSL can successfully calculate the transport and distribution of humidity in the RSM.One possible reason for why the improvement of the rainfall simulation by the NDSL scheme is not very significant is that the selected case in this study is not a heavy rainfall case.For a heavy rainfall case, the large discontinuity of a humidity field is expected, which means higher possibility of negative value occurrences in the original IsoRSM.Further studies will continue to examine how the NDSL can improve skills for the precipitation simulation in a heavy rainfall cases.

Summary and discussion
In previous studies, the RSM was utilized to simulate the Fukushima Daiichi nuclear power plant accident.The results exhibit severe noise in the simulated radioactive tracer fields (i.e., iodine-131 and cesium-137).This noise is due to the Gibbs phenomenon, wherein discontinuities in positivedefinite fields create negative values after spectral transformations.This problem is common in spectral model systems.The spectral tracer advection is replaced with a semi-Lagrangian advection scheme to prevent the Gibbs phenomenon in the tracer output.Prognostic tracer fields are  calculated using the semi-Lagrangian method and are only considered in grid space.The Gibbs phenomenon does not occur when the spectral space transformation is not performed.
The semi-Lagrangian method used in this study is the NDSL scheme, which has the advantages of efficiency and simplicity.Because the NDSL has been previously applied in a global model system only, a regional version of the NDSL is developed in this study.For this application, the boundary conditions are applied by defining simple weighting functions.The regional version of the NDSL scheme is verified by performing idealized experiments using horizontal and vertical advection.These idealized experiments transport a particular disturbance to the uniform wind field.The results show that the shape is well maintained and the mass conservation is satisfied during advection.Therefore, the regional version of the NDSL is successfully applied in the RSM.
Two experiments are performed for the Fukushima case study to evaluate the NDSL advection scheme in the RSM.The ORG experiment is performed by the original RSM, and the SL experiment is produced by the NDSL version of the RSM.The ORG run shows severe errors in tracer fields induced by the Gibbs phenomenon.Errors appear as a ringing signal that extends zonally and meridionally from the emission point.Additionally, relatively strong ring-shaped noise is captured around the emission point.This noise is clearly removed when the tracer advection component is replaced by the NDSL scheme.The SL run shows that the NDSL advection scheme can capture the major transport of tracers without any noise from the Gibbs phenomenon.This finding is the clearest advantage of the NDSL scheme in the tracer field simulation.In the humidity field, the ORG experiment produces some negative values in the lower troposphere.However, the SL experiment does not exhibit such negatives; both SL and ORG capture the detailed distribution of the humidity field.The precipitation fields from the ORG and SL ex-periments are similar, which means that the NDSL properly calculates the humidity field.
This study reveals that replacing the tracer advection scheme with a semi-Lagrangian scheme can eliminate the Gibbs phenomenon in a regional spectral model.However, the simulated surface depositions of radioactive tracers still deviate from the observations, and precipitation from the SL experiment does not show significant improvement, even though the NDSL removes severe errors.Note that the objective of this study is to determine the feasibility of the NDSL advection scheme in a regional spectral model.Thus, some quantitative validations from experiments are not included in this study.These results may be improved upon by applying enhanced physical parameterizations and advanced formulas for tracer surface deposition processes.

Figure 1 .
Figure 1.Schematic of two-dimensional advection for (a) the traditional backward semi-Lagrangian scheme and (b) the NDSL scheme."D," "M," and "A" indicate the departure point, midpoint, and arrival point, respectively.The bold arrow indicates the wind vector at the midpoint.

Figure 2 .
Figure 2. (a) The structure of the lower-left corner in the regional domain for boundary treatment.(b) Weighting function of the global base field (black) and regional model field (gray) along the x axis, which includes the inner domain.This example assumes that the grid size of the boundary and buffer zones is 3.

Figure 3 .
Figure 3. Result from the idealized horizontal advection experiment.(a) Virtual concentration every 1000 time steps by the NDSL advection scheme and (b) results at the 3000th time step from the NDSL (red contour) and the analytic solution (black contour).The shaded values indicate the differences between the NDSL results and the analytic solution.

Figure 4 .Figure 5 .
Figure 4. Virtual concentration from the idealized vertical advection experiment.The red and black lines indicate the transported concentration by the NDSL and the analytic solution, respectively.