The impact of periodization methods on the kinetic energy spectra for limited-area numerical weather prediction models

The paper deals with the comparison of the most common periodization methods used to obtain spectral fields of limited-area models for numerical weather prediction. The focus is on the impact that the methods have on the spectra of the fields, which are used for verification and tuning of the models. A simplified model is applied with random fields that obey a known kinetic energy spectrum. The periodization methods under consideration are detrending, the discrete cosine transform and the application of an extension zone. For the extension zone, three versions are applied: the Boyd method, the ALADIN method and the HIRLAM method. The results show that detrending and the discrete cosine transform have little impact on the spectra, as does the Boyd method for extension zone. For the ALADIN and HIRLAM methods, the impact depends on the width of the extension zone – the wider the zone, the more artificial energy and the larger impact on the spectra. The width of the extension zone correlates to the modifications in the shape of the spectra as well as to the amplitudes of the additional energy in the spectra.


Introduction
As the horizontal resolution of numerical weather prediction (NWP) models steadily increases, rigorous evaluation of the models' kinetic energy spectra at small scales becomes more important.At horizontal scales between a few hundred and a few km, which are the focus of mesoscale NWP models, statistical properties of turbulent motions in the upper troposphere and lower stratosphere have been shown to be consistent with the Kolmogorov scaling of turbulence defined by the energy dissipation rate.Scale distribution of turbulent motion energy is described by the k −5/3 law, where k is the horizontal wave number.The scaling is horizontally nearly isotropic as it was confirmed in aircraft observations (e.g.Nastrom and Gage, 1985;Lindborg, 1999) and numerical model outputs.In models, kinetic energy spectra at the shortest scales, represented by the model grid resolution, result from a complex interplay between the parameterizations and the numerical schemes of the models.In particular, horizontal diffusion schemes strongly affect the shape of the spectra at small scales and the models' ability to simulate turbulent processes (e.g.Skamarock, 2004).As a result, the effective model resolution is on average several times coarser than defined by the model numerical truncation (e.g.Frehlich and Sharman, 2008).For example, for the ALADIN model applied for operational NWP in a number of countries across Europe, the effective resolution was found to be about 6 x, where x is the grid distance of the model (Blažica et al., 2013).
Comparison between the kinetic energy spectrum in NWP models and the expected k −5/3 scaling is a way to assess how well small scales are resolved by the model.Frehlich and Sharman (2008) evaluated spectra for several NWP models with horizontal resolutions of about 10 km at levels close to the aircraft data and found that spectra start to deviate from the observed slope already at several hundred km.The slope of the spectra is significantly different for the rotational and divergent components of kinetic energy with the latter more shallow at all scales (Blažica et al., 2013).
In most studies, calculation of kinetic energy spectra is based on the Fourier transformation.An alternative, the computation of structure functions (Frehlich and Sharman, 2004) as performed in Frehlich and Sharman (2008), is typically not applied by the NWP community.In particular, for the large community of users of spectral limited-area models in Europe, which are based on the application of the Fourier transform, the use of readily available spectral coefficients of the model prognostic variables has practical advantages.
In Europe, two spectral limited-area models (SLAMs) used for operational NWP are the Aire Limitée Adaptation Dynamique Développement International (ALADIN, http: //www.cnrm.meteo.fr/aladin/)and the spectral version of the High-resolution limited-area model (HIRLAM, http://www.hirlam.org).The advantage of SLAMs compared to gridpoint models is the accuracy and computational efficiency of the high-order derivatives (Termonia et al., 2012).The application of the Fourier transforms in these models requires bi-periodicity of the prognostic fields which is in both models achieved by the so-called extension zone or E zone (Haugen and Machenhauer, 1993).The application of the E zone adds an additional belt of grid points in both x and y directions, in which the fields are periodized using trigonometric (HIRLAM) or spline functions (ALADIN).Recently, the Boyd method (Boyd, 2005) was employed to the AL-ADIN model as well (Termonia et al., 2012;Degrauwe et al., 2012); the method uses infinitely differentiable windowing functions in the extension zone.
Two other periodization methods, which do not require the extension zone, are detrending (Errico, 1985) and the discrete cosine transform (DCT, Ahmed et al., 1974;Denis et al., 2002).The detrending method removes large-scale trends in each line and in each column of the grid-point field and it has been widely used for the computation of spectra.The DCT method creates a symmetric field by taking a mirror image of the original function.Neither of these two methods is suitable for the forecasting as they alter the values or the dimension of the prognostic variables.Instead, they are used a posteriori for the verification and tuning purposes (e.g.Skamarock, 2004;Ricard et al., 2013;Brousseau et al., 2012).
In the case of the ALADIN model and its versions AROME (Seity et al., 2011) and ALARO (Gerard et al., 2009;De Troch et al., 2013), the presence of the E zone in the spectral computation assumes that the kinetic energy, added within the extension zone, does not have a significant impact on presented spectra.While this may apply to some levels and situations, it does not apply in general, as seen in Fig. 1 which compares kinetic energy spectra in the stratosphere in two experiments with everything the same except the width of the E zone.It can be seen that the energy spectrum for the E zone with a width around 50 km (i.e.11 grid points with x = y = 4.4 km) is characterized by an energy increase at the scale corresponding to the width of E zone.A wider E zone which includes 25 points (i.e.width around 100 km) shows increased energy starting around 100 km.The difference between the impact of the two E zone widths disappears above about 250-300 km and under 50 km.Figure 1 also illustrates the impact of the selected method, i.e. application of the E zone or detrending, on the large-scale portion of the spectra, demonstrating that the two spectra cannot be compared quantitatively.Further details on this aspect are available in Blažica et al. (2013).
The present study deals with the impact of various periodization methods on the shape of kinetic energy spectra over a limited area.The focus is on spectral limitedarea modelling as applied in the NWP models ALADIN and HIRLAM.We present a systematic evaluation of several methods of periodization which are in use as a function of the domain size and the width of the extension zone.In order to allow a clean comparison, we use a simple model and an idealized flow with a known spectrum.
Details of the testing methodology and selected periodization methods are presented in Sects. 2 and 3, respectively.Section 4 discusses the resulting spectra while conclusions are stated in Sect. 5.

Overview of the periodization techniques
The following periodization methods are taken into consideration: the extension zone with trigonometric functions (as used in the spectral version of the HIRLAM model), the extension zone with spline functions (as used in the ALADIN model), the extension zone with the Boyd method (a newer option for the ALADIN model), the detrending method and the discrete cosine transform.The first three methods keep the basic field unchanged and perform periodization in the artificial, extended zone, which makes them appropriate for use in spectral modelling.The latter two methods are only appropriate in post-processing for spectral analysis of the model outputs.

The extension zone in the HIRLAM and ALADIN models
The extension zone (Haugen and Machenhauer, 1993) is essential for the spectral models HIRLAM and ALADIN in order to make use of fast Fourier transforms (FFTs) possible.The grid-point space calculations are carried out over a domain consisting of N xi × N yi horizontal grid points.In order to obtain bi-periodic variations needed for Fourier transforms and spectral space calculations this domain is extended to N x × N y horizontal grid points (Fig. 2).
The current application of ALADIN is based on 11point E zones at the northern and eastern boundaries while HIRLAM uses a more flexible extension depending on the actual modelling or assimilation problem.The biperiodization technique to obtain the values of prognostic fields in the extension zone is only applied for calculation of initial and lateral boundary conditions for the spectral model.During the data assimilation, no bi-periodization is needed since the extension zone is treated like other areas without observations.
The extension of the prognostic fields is performed by using spline (ALADIN) or trigonometric (HIRLAM) functions, which are applied first row by row and second column by column.
The formulas for the row-by-row part in ALADIN are with M and Here, W is a prognostic variable, N is the width of the extension zone (N = N x − N xi + 1) and N i is the index of the current point within the extension zone.
For HIRLAM, the formula is where N and the coefficients g j,k (k = 0, 1, 2 and 3) are determined to make the extrapolation fit the original grid-point values exactly for i = 1, 2, N xi − 1 and N xi .
In ALADIN, solutions over the extension zone are additionally smoothed by a two-dimensional 9-point weighted average for each point in the E zone.One motivation for the smoothing was the removal of small-scale dynamical features in the interior domain that are reflections of small-scale features in the extension zone caused by the isotropic truncation.
Because the E zone is considered to have a negligible impact on model calculations, it is usually included in the kinetic energy spectra of dynamical fields and in the spectra of the forecast-error variances used in data assimilation.As discussed in Blažica et al. (2013), the kinetic energy spectra with the extension zone included are regularly used in AL-ADIN and HIRLAM studies of the model performance (e.g.The impact of periodization on the energy spectra Horvath et al., 2011) and for tuning of the model numerics (e.g.Váňa et al., 2008) and physics (e.g.Bengtsson et al., 2012).
In the HIRLAM model, imbalances in the extension zone as a possible source of problems for the initialization are treated by geostrophic balancing of the gravest vertical modes (tendencies of the fastest gravity modes are put to zero).In HIRLAM, the extension is carried out for the lateral boundary data only.

The extension zone based on the Boyd method
The Boyd method (Boyd, 2005) is another method that uses the extension zone, but instead of applying a spline or a trigonometric periodization operator, the method makes use of the prognostic variable values from outside the model domain.This is not seen as a problem because periodization is done on the data from the host model.The Boyd method can be seen as a decrease in size of the extension zone, but the transition from the solution on one side of the domain to the solution on the other side still introduces some aliasing.
The values in the extended area are a combination of the additional values beyond the eastern/northern boundary and the additional values beyond the western/southern boundary.The transition between the two is achieved by a windowing function (Boyd, 2005;Termonia et al., 2012): Here the interval to be periodized is [− , ], the physical domain is limited to the smaller interval [−χ , χ ], L is a tunable parameter (in our experiment set to 1.6) and erf is the error function.The variable is first multiplied with the windowing function and next a summation of the shifted products is performed: This method has recently been adapted to the ALADIN code and can now be run operationally.Its performance in the HARMONIE model was demonstrated by Termonia et al. (2012) and Degrauwe et al. (2012), who showed that the Boyd method outperforms the spline-based periodization method, especially in the cases of strong dynamical forcing at the lateral boundaries.This suggests that the method applied to make fields periodic has an impact on the model performance, for example through the spectral calculation of spatial derivatives during the model integration.

The detrending method
The detrending method (Errico, 1985(Errico, , 1987) ) removes the scales that are larger than represented by the limited domain and thus not resolved.It consists of the computation of the linear trend between the end points along each row and column of data, followed by the removal of the trend: Here, s j denotes the linear trend and W is the modified field.The procedure matches the first and the last point of the row/column, thus making the fields periodic.The manipulated field usually contains a pattern of lines, created by the linear treatment.As mentioned by the author, the method should not be used if the fields are noisier at the boundaries than in the interior of the domain.

The discrete cosine transform
The discrete cosine transform (Ahmed et al., 1974;Denis et al., 2002) is a method most commonly used for image processing.The periodization is achieved by adding a mirror image to the variable field.Due to the even symmetry of the created field, the Fourier transform consists only of the cosine functions, hence the name of the method.
The periodic function is constructed by taking the position j = −1/2 as a mirror: a discrete Fourier transform is then applied to the new function, centred on i = −1/2.While this transform method is widely used for compression of digital images, its use for atmospheric spectral analysis has only recently received attention.It has recently been applied for the computation of kinetic energy spectra in AL-ADIN as it avoids problems associated with the extension zone (Ricard et al., 2013;Brousseau et al., 2012).

Demonstration of the above methods
For a demonstration of the periodization procedure, all the methods were applied to an arbitrary field with a small domain and a large extension zone (60 × 60 grid points, with 18 points in the extension zone).The detrending method, which does not require an extension zone, was applied on the entire domain (60 × 60 grid points).Figure 3 demonstrates periodization of one grid-point line for each of the discussed methods (except DCT).The two-dimensional prog- nostic fields after the application of the periodization methods are presented in Figs. 4 and 5.
The properties of the fields with the Boyd and detrending methods for periodization look similar to the original field (top left in Fig. 4) while the HIRLAM and ALADIN methods produce patterns along rows and columns of the E zone (which are partially omitted after the ALADIN smoothing procedure) and an area with absolute maximum in the top right corner of the domain (where the values are obtained after extrapolation in both horizontal directions).Both phenomena impact the spectra, as will be shown below.

Simulation of random wind fields obeying a known kinetic energy spectrum
The experiment to test the impact of the methods presented above consists of applying these algorithms to random fields, generated in a way that they correspond to wind fields obeying a known kinetic energy spectrum.Such random fields are constructed as follows.We may represent the wind field in spectral space by the Fourier series coefficients ûkl , vkl where k and l are the horizontal wave numbers in the x and y direction, respectively.It is then possible to generate a random wind field obeying a specified kinetic energy spectrum S(k * ) where k * = √ k 2 + l 2 is the onedimensional wave number.Such a random wind field is, for example, given by where Parameters G Re u , G Im u , G Re v and G Im v are normally distributed (Gaussian) random numbers with mean = 0 and standard deviation = 1.The factor 2π k * follows from an integration over all wave numbers obeying k * = √ k 2 + l 2 to obtain the total kinetic energy for this particular one-dimensional wave number.
For the presented experiment the given spectrum was set to S(k * ) = (k * ) −5/3 and the domain size (N x , N y ) to 432 × 432 grid points.The created spectral fields were transformed to grid-point space by a two-dimensional inverse FFT, where for the cases with the extension zone the data matrix is reduced by the width of the extension zone (to N xi , N yi -see Fig. 2).Note that the outer domain remains the same for all the discussed methods; thus the wider extension zone means a reduction of the inner, physical field.After the application of the selected method for periodization, the fields are transformed back to spectral space by a two-dimensional direct FFT and their new kinetic energy spectrum is compared to the original one.
The applied simulation method favours the detrending method as artificially constructed fields contain nearly no trend to be removed.In reality, fields can have a large trend (e.g.north-to-south temperature gradient, wind gradient when a trough enters the domain etc.) and detrending could produce a much larger effect.Mind also that for the Boyd method the grid-point values are obtained from the "true" spectrum with the correct theoretical scaling in the extension zone.This is not the case in the lateral-boundary couplings of real models where the goal is to nudge the solution to the one of the host model grid-point values as best as one can, obeying the spectrum of the host model.

The impact of periodization on the kinetic energy spectrum
In this section we discuss the kinetic energy spectra for all the above mentioned methods.The methods are named after the models that use them.The results, based on averaging over 500 realizations of the random wind field, are depicted in Fig. 6.For the methods that make use of the extension zone, the spectra for two different widths of the extension zone (18 and 48 grid points) are presented.For the ALADIN model, we show two sets of spectra -one after the spline function periodization (referred to as ALADIN) and the other after applying the smoothing to the periodized fields (as done operationally; referred to as ALADIN SMOOTH).
The resulting spectra show significant differences among the methods under investigation.While the Boyd and the detrending methods remain very close to the original spectrum, there is a significant increase of energy for the ALADIN and   HIRLAM methods.The similar shapes of spectra for these two methods show that the deformation is mainly a consequence of the extension zone geometry and not of the function used for the periodization.For the extension zone width of 18 grid points the added (excess) energy increases gradually towards shorter scales.For the 48-point zone, this process is also observed, but the most distinctive feature of the spectra is a bump, centred approximately at wave number 8 with an amplitude about one order of magnitude larger than for the original spectrum.A more careful inspection reveals that a very small bump can already be seen in the 18-point E zone for the HIRLAM method, centred at about wave number 25.We investigate the location, amplitude and cause of this phenomenon later in the paper.
In the wide extension zone, the ALADIN SMOOTH spectrum first follows the ALADIN spectrum in the bump region, but descends from it and approaches the original spectrum at about wave number 50.Towards the shortest scales (at about wave number 60), the amplitude of the spectrum becomes even smaller than the original k * −5/3 because the smoothing operator filters out some of the short-scale waves.The slope of the spectrum for the ALADIN SMOOTH method is thus steeper than the original one.The curve resembles Fig. 1, where the ALADIN SMOOTH method was applied to real forecast fields.
Figure 7 shows the ratio between the computed spectra and the prescribed k * −5/3 spectrum for the detrending method, the DCT, and for several widths of the extension zone for the Boyd method.At large scales, all spectra contain energy which exceeds the energy of the original field.All methods provide spectra that slowly approach the prescribed spectrum as scales become shorter.The energy ratio remains stable for wave numbers above 30.The spectrum of the detrending method remains above the amplitude of k * −5/3 spectrum at all scales while the Boyd method loses some of the initial energy, depending on the size of the extension zone.The best match to the original spectra presents the DCT method.The differences among the spectra are small compared to the impact of the ALADIN and HIRLAM methods, on which we focus next.
For the HIRLAM and ALADIN methods, the spectra with several widths of the extension zone are shown in Fig. 8.The deformation of the spectra increases and moves upscale with increasing size of the E zone, supporting the use of a narrow extension zone for a smaller impact on the spectra.The location of the bump shifts nearly linearly towards smaller wave numbers.This confirms that the location of the bump is closely related to the width of the extension zone.The AL-ADIN SMOOTH spectra resemble the ALADIN spectra at large scales, but the energy decreases significantly for shorter scales due to the smoothing effect.
A similar ratio as in Fig. 7 is shown in Fig. 9 for the HIRLAM and ALADIN spectra.Because for these two methods there are no special requests for the size of the extension zone, periodization was done for 25 widths, from 0 to 216 grid points, increasing the width of the extension zone by 8 grid points.The size of the extension zone in percentages of the full field is shown on y axis.While such large sizes of the extension zone are not used in practice, these plots emphasize the importance of the extension zone width to the location of the spectral bump and to an overall deformation of the spectrum.
To further investigate the large energy values in the AL-ADIN and HIRLAM spectra, we perform row-wise onedimensional FFTs in the x direction.Thus obtained spectra are averaged over rows 1 to N yi (384) separately and over rows N yi + 1 (385) to N y (432) separately (see Fig. 10).The former spectra mainly contain the information in physical space with only the last 48 points of each row in the extension zone and should thus enclose the effect of periodization in single rows.The latter spectra on the other hand only contain the information in the extension zone and should enclose the effect of performing periodization line-by-line.Additionally, these spectra also include large values of the extended fields in the top right corner (see Fig. 4).
Averaging the one-dimensional spectra in the physical field (with only the last 48 grid points reaching into the extension zone) shows that the added energy only occurs at scales larger than the width of the extension zone (Fig. 11).The reason for this is that the spline or trigonometric functions complete the series of data with shapes varying from a single wave towards a flat line, which therefore projects the added energy on the scales which correspond to waves with scales larger than the extension zone width.Smoothing the ALADIN method fields further reduces the added energy.
Averaging over the rows in the extension zone gives much larger amplitudes of the spectra across all scales.The contribution from large scales to the spectra is larger than in the previous case and it is not affected by the smoothing operator.It is most likely that the cause is the top right corner of the extended domain, which was not included in the physicalfield spectra.There is also some additional energy towards shorter scales, which was not present in averaging in physical space.This small-scale part of the spectra most likely comes from the line-by-line periodization (extrapolation), the effect of this was not present in the physical-field spectra.This shorter-scale excess energy is however efficiently removed by applying the smoothing filter in ALADIN SMOOTH.
The analysis of one-dimensional spectra obtained by the ALADIN and HIRLAM methods enables a valuable insight into the structure of the spectra.We showed that the bump in the spectra at large scales occurs for two reasons: one is the periodization function, which completes the data with shapes that vary from one wave to a flat line; the other is the extreme values in the top right corner of the periodized domain, which only influence the spectra in the extension zone.The shortscale excess energy comes from row-by-row and column-bycolumn periodization of the variables.The smoothing operator efficiently filters out the short-scale noise and some of the large-scale energy.
The geometry of the domain in this experiment is very similar to that of ALADIN model which provided data used to construct Fig. 1.Comparison of these simple spectra to the spectra of real forecast fields where the ALADIN SMOOTH method was applied helps to interpret the results with the real fields: the smoothing operator efficiently removes the excess short-scale energy, added by the extension zone, and for this reason the spectra of the extension zone method and the spectra of the detrending method match below the scale of 50 km.The large-scale bump, originating from the treatment of the top right corner of the domain, however still remains, and causes the deviation of the extension zone spectra from the detrended spectra.

Conclusions
We presented the kinetic energy spectra of a limited area, obtained with different periodization methods.In particular, we focused on the methods that make use of the extension zone, as in the ALADIN and HIRLAM models, and we investigated the impact of the extension-zone widths on the structure of the spectra.A narrow extension zone is desired for a lower computational cost while the application of larger extension zones is preferable in variational data assimilation since horizontal correlations of background errors are computed in the spectral space including the extension zone.In the case of too narrow extension zones there may be an artificial influence of observations situated close to one of the lateral boundaries on grid points close to the lateral boundaries on the opposite side of the domain (Gustafsson et al., 2001).It is thus better to use a wider extension zone for data assimilation and to switch to a narrower one for bi-periodization of prognostic fields and model integration.
The experiments presented in this study offer the following conclusions: -The spectra obtained by the discrete cosine transform, the Boyd and the detrending method match the simulated idealized spectra well.The Boyd method is from this point of view therefore a suitable method to replace the existing spline interpolation in the ALADIN model, in agreement with the study of Degrauwe et al. (2012).
The detrending approach and the discrete cosine transform are more suitable for the periodization for diagnostic purposes, because, contrary to Boyd method, no additional data from outside the model domain is required.
-While the impact of the HIRLAM and ALADIN methods on energy spectra is small in cases with a narrow ex-  tension zone, it increases as the E zone becomes wider, deforming the shape and thus also the slope of the spectra.The energy spectrum including the E zone information is characterized by a bump, the location of which is closely related to the width of the extension zone.The contamination of the spectrum by the E zone at small scales is a consequence of treating the data matrix lineby-line and it is efficiently removed by the smoothing filter, although even after smoothing the spectrum slope is not correct.
-In many HIRLAM and ALADIN applications, the impact of periodization may pass nearly unnoticed in the spectra as extension zones are usually narrow and prognostic fields have large amplitudes at scales most affected by the E zone.This primarily applies to the fields in the planetary boundary layer and to a smaller extent to fields in the troposphere.It does not apply to stratospheric circulation dominated by large-scale waves, as illustrated in Fig. 1.
-The original HIRLAM and ALADIN algorithms should not be applied for calculation of spectra with wide extension zones since the constraint on preservation of normal gradients together with the selected extrapolation functions leads to strong amplification of variations in the extension zone.

Figure 1 .
Figure 1.Kinetic energy spectra of the ALADIN/SI model with 4.4 km horizontal resolution, averaged over model levels between 5 and 250 hPa.Top: comparison of ALADIN spectra with 11 and 25 points in the E zone.Bottom: comparison of ALADIN detrended spectra and spectra with E zone.The dotted lines have slopes (k * ) −3 and (k * ) −5/3 where k * is the one-dimensional wave number.

Figure 2 .
Figure 2. A typical setup of the NWP model domain with the extension zone.

Figure 3 .
Figure 3. One-dimensional periodization of an arbitrary wind field for the discussed methods except for DCT.The thin grey lines denote the borders of the extension zone.

Figure 4 .
Figure 4.The arbitrary 2-D wind fields before (top left) and after the periodization for the HIRLAM, ALADIN, ALADIN SMOOTH, Boyd and detrending methods.The black lines denote the borders of the extension zone.

Figure 5 .
Figure 5. Same as in Fig. 4, but for DCT.The periodization is shown only for x direction.

Figure 6 .
Figure 6.Kinetic energy spectra for all the discussed techniques.The results with two widths of the extension zone are presented for HIRLAM, ALADIN and Boyd methods: a narrower one with 18 points (top) and a wider one with 48 points (bottom).

Figure 7 .
Figure 7.The ratio between the computed spectrum and the prescribed k * −5/3 spectrum for the detrending method, the discrete cosine transform and the Boyd method for several widths of the extension zone.

Figure 8 .
Figure 8. Kinetic energy spectra for HIRLAM (top), ALADIN (middle) and ALADIN SMOOTH (bottom) methods for several widths of the extension zone.

Figure 9 .
Figure 9.The ratio between the periodized spectrum and the original k −5/3 spectrum for several widths of the extension zone for HIRLAM (top), ALADIN (middle) and ALADIN SMOOTH (bottom) methods.The ratio between the size of the extension zone and the full domain is shown on y axis.

Figure 10 .
Figure 10.An arbitrary grid-point field with N xi = N yi = 432 and the extension zone width of 48 points, periodized with the HIRLAM method.

Figure 11 .
Figure 11.Average one-dimensional spectra only in the extension zone (EZONE, solid lines) and only in the physical field (PF, dashed lines) for HIRLAM, ALADIN and ALADIN SMOOTH methods.