## Introduction

Hyper-viscosity is widely used in numerical models for the weather and climate prediction as a viscosity (or diffusion) scheme (Jakob-Chien et al., 1995; Gelb and Gleesen, 2001; Cheong et al., 2004; Krishnamurti et al., 2006). Most commonly used hyper-viscosity formula is represented with the fourth order differential operator of Laplacian-type (e.g., Jakob-Chien et al., 1995). The Weather and Research Forecasting (WRF) model (Skamarock et al., 2008), one of the most popular community model for the atmospheric sciences, adopts a sixth-order hyper viscosity or the third-order Laplacian. As the order of the viscosity becomes higher, a sharper cutoff (*or* filtering) of the small-scale disturbances can be achieved while leaving the largescales affected little (Haltiner and Williams, 1980; Gelb and Gleesen, 2001; Cheong et al., 2004). Discretization of the high-order viscosity with an explicit time scheme as in the WRF model leaves some undesirable properties such as restriction to the viscosity coefficient and a large number of stencil points: The number of stencil points in one coordinate direction is given by 2*q*-1 with *q* being the order of hyper-viscosity.

Implicit treatment of the hyper-viscosity is free from the restriction to the filter coefficient and can be used as a high-order spatial filter (Park et al., 2011; Cheong and Jeong, 2015; Kang and Cheong, 2017) as well. The high-order filter is known to provide the same performance as multiple (typically more than 30 times) successive operation of the explicit hyper-viscosity (Lauritzen et al., 2015; Kang and Cheong, 2017). The high-order filter can be implemented by various discretization methods provided that the Helmholtz equation, which constitutes the building block of the high-order filter, can be discretized. The high-order filters in the previous studies are constructed by incorporating the Fourier spectral method in the zonal direction and the Galerkin method or finite difference method (FDM) in the meridional direction (Park et al., 2011; Cheong and Jeong, 2015; Cheong et al., 2015). One disadvantage of using the high-order filters in the previous studies may be the inefficiency due to the Fourier transform needed for the longitudinal discretization. Although the Fourier spectral method for those filters based on the Fourier-finite element method (FFEM) as in Heinrich (1996) and Cheong et al. (2015) provides an accurate discretization as well as an exponential error-convergence rate (Canuto et al., 1988), the overall accuracy is determined by the finite element method (FEM) used for the meridional direction. In this regard, the high-order filter with two dimensional FEM may be more useful than those in the previous studies, because it does not need the Fourier transform or the Fast Fourier Transform (FFT) module which is not considered well suited for the parallel computing on distributed memory (e.g., Cortese and Balachandar, 1995). Another important advantage over the previous filters may be that the FEM filter can be implemented on arbitrary grids.

In this study, the high-order filter with quadrilateralbasis FEM on the spherical-surface limited area domain is developed for a use in the WRF model as a hyperviscosity as well as a filter to separate a given field into two parts having different horizontal scales. This paper is organized as follows. In the following section, the procedure to build filter matrices using the FEM with quadrilateral basis functions is described. Section 3 provides the evaluation of the high-order filter in terms of the normalized error and the error convergence rate. The numerical simulations of the tropical cyclones using the WRF model with the high-order filter as a hyper-viscosity are given in section 4. The final section presents the summary and conclusions.

## Two Dimensional Finite Element Method for the High-order Filter on the Spherical-surface Local Domain

The spherical Laplacian operator on a limited area domain can be written as (Haltiner and Williams, 1980; Jakob-Chien et al., 1995; Krishnamurti et al., 2001):

where *λ* and *θ* represent the longitude and the latitude, respectively, and *λ*_{1} and *λ*_{2} (*θ*_{1} and *θ*_{2}) are the westand east- (south- and north-) boundaries of the domain, respectively (Fig. 1). Helmholtz equation, which is the first-order Laplacian-type filter, is written as (e.g., Cheong et al., 2015)

where *ψ* and *F* are scalar functions which are square integrable on the spherical local domain. The finite element procedure will be explained below to solve Eq. (2) for a given forcing function *F*, where the boundary condition is imposed as a gradient-vanishing condition (Denis et al., 2002; Park et al., 2011).

An arbitrary function on the sphere is represented as a finite sum of basis functions defined on a local domain as (Haltiner and Williams, 1980; Durran, 1999)

where *ψ _{i,j}* means grid-point data,

*φ*implies the basis function, and

_{i,j}*K*and

*N*are the total number of grids in the longitudinal and latitudinal direction, respectively. Each basis function consists of four shape functions given as

where *λ _{i}* and

*x*are the grid points in the longitude and latitude, respectively, Δ

_{i}*λ*[=(

*λ*

_{1}–

*λ*

_{2})/

*K*] is the grid interval in the longitude, and = (Δ

*x*)

_{j}=

*x*–

_{j}*x*= sin

_{j-1}*θ*– sin

_{j}*θ*. The grid interval in the meridian is set a constant value, i.e., Δ

_{j-1}*θ*= (

*θ*

_{2}–

*θ*

_{1})/

*N*. The basis function of

*φ*

_{i,j}is illustrated in Fig. 2. The shape function has a local peak at the grid point

*φ*

_{i,j}and decreases monotonically toward the adjacent grid boundaries. Note that, unlike the one-dimensional shape function (e.g., Durran, 1999; Cheong et al., 2015), the function shows a quadratic behavior except on the edge of grid-boxes. It is clear that the function

*φ*

_{i,j}is associated with adjacent 8 grid points, yielding 9 point stencil for the FEM discrete analogue of Eq. (2).

The basis functions are multiplied to the elliptic equation (2), and integrated in the model domain (Haltiner and Williams, 1980; Heinrich, 1996; Durran, 1999; Cheong et al., 2015), e.g.,

Applying the integration by part to Eq. (5), a matrix equation is obtained.

where **D** and **F** are sparse matrices, with only 9 nonzero entries, having the size of *KN×KN*, and * ψ* and

**F**are column vectors with the length of

*KN*consisting of the grid-point values of

*ψ*and F, respectively. The matrices are diagonally dominant irrespective of the magnitude of the coefficient

*ν*, therefore inversion of

**D**can be done without singularity. As demonstrated in Kang and Cheong (2017), Eq. (6) can be solved in two ways. Of the two methods, the inverse matrix method is more efficient than the row reduction method, because it requires only forward operation. In this study, we will not pay attention to the computational efficiency in association with the high-order filter. This will be pursued in depth in a subsequent paper together with a parallel implementation on a distributed memory high-performance computers. At the moment, Eq. (2) is solved by the row reduction method with a focus on the accuracy and the applicability to the WRF model.

The high-order filter equation is the same as in the literature (Cheong et al., 2004; Cheong and Jeong, 2015; Kang and Cheong, 2017): The spherical highorder filter of order *q*, for the given function *ψ*, is written as

where * means the function obtained by filtering. As shown in Cheong et al. (2004, 2015), the filter coefficient in the time-integration of models is determined in such a way that the smallest scale is dampened by half at a prescribed time scale. On the other hand, in the case of scale decomposition for an arbitrary grid-point data on the spherical surface such as the separation of the disturbance scale and the large-scale environmental field, the filter coefficient should be specified to damp down the amplitude of the target scale (*or* filter scale) by half. To see how the matrix size is varied with the filter order *q*, one may refer to the recent study of Kang and Cheong (2017).

## Evaluation of the Accuracy

The spherical harmonic function of *P*^{8}_{9}cos8*λ* is used as a reference field to evaluate the accuracy of the filter, where *P*^{8}_{9} is the associated Legendre function of the degree 9 and order (*or* zonal wavenumber) 8. The limited area domain is set as *λ*_{1}=0°E, *λ*_{2}=45°E *θ*_{1}=19.47°S and *θ*_{2}=19.47°N, and filtering is performed for 16 cases with different filter order and the filter coefficient. The filter coefficient is given as *ν* =10^{−4} or *ν* =10^{−2}, which can damp down the amplitude of the forcing function (*P*^{8}_{9}) by (1+*ν*90^{q})^{−1}. The number of gridpoint is set as *K*=*N*=20×2^{s} with s an integer between 0 and 3. The accuracy is evaluated in terms of the root mean squared error (RMSE) and the error convergence rate (Fig. 3) for the different filter order (left panel of Fig. 3) and for the different filter coefficient (right panel of Fig. 3). The RMSE appears to vary in the range of *O*(10^{−2})~*O*(10^{−4}), and the convergence rate is the second-order, as is expected. The convergence rate is the same as the global FFEM elliptic solver (Cheong et al., 2015). The errors decrease with both the decreasing filter order and viscosity, and they appear to be the same order as those in the FFEM method. The fact that the accuracy as well as the convergence rate is not affected by the change of discretization method in the zonal direction reflects that the overall numerical performance is determined mainly by the discretization method in the meridional direction. As consistent with the previous study with the FFEM, the errors are sensitive to the order of the filter as can be seen in the left panel of Fig. 3.

Figure 4 provides two examples of filtering used for decomposing a limited-area scalar field into the largescale and small-scale components. The scalar field is composed of the zonal-mean component and two nonzonal disturbances of Gaussian-bell structure. In these examples, the reference scale (*or* filter scale) is specified as *M _{c}*=6 or 24 of the spherical wavenumber (i.e., the degree of spherical harmonic function), and the filter order is either 1 or 4. It turns out that the filtered fields have quite different structure one another depending on these two parameters. These results suggest that separation of disturbances from the given field may be done successfully as intended only when the filter scale (or the scale of disturbances which should be removed by filtering) is given appropriately. For instance, if one wants to separate the tropical cyclones from the observed meteorological field in the bogussing procedure (e.g., Kurihara et al., 1993; Kwon and Cheong, 2010), the information on the size of the tropical cyclone under consideration is essential. For this purpose, Kwon and Cheong (2010) and Cheong et al. (2011) used two parameters, the radius of 30kt and the maximum wind speed, and the filter scale was given to be between 14 and 18.

## Application to the Tropical Cyclone Prediction

Tropical cyclone prediction using numerical model requires a sophisticated initialization technique in order to represent appropriately the structure of the tropical cyclone (Kurihara et al., 1993; Kwon and Cheong, 2010; Yoshiaki et al., 2014). In this study, the bogus method developed by Kwon and Cheong (2010) is used, in which the tropical cyclone is represented with the 3-dimensional empirical formula. The observed tropical cyclone is known to be well approximated by gradient wind balance (Tuleya and Kurihara, 1974; Kurihara et al., 1993; Persing et al., 2013; Leonov, 2014; Yoshiaki et al., 2014), and hence this constitutes a key factor in the empirical formula of the bogus typhoon. The primary goal of the bogussing, proposed originally by Kurihara et al. (1993), is to replace the coarse-resolution tropical cyclone in the analysis data by a high-resolution realistic vortex. The high-order filter explained above is incorporated in the simulations with two different purposes: One is to separate the disturbance of tropical-cyclone scale from the largescale environmental field in the initialization procedure, and the other is to use as a hyper-viscosity to get rid of grid-scale numerical noises generated by the nonlinear terms. The WRF v3.4.1 model is used for numerical simulations with numerical options which are the same as in Kang et al. (2017). The model is time integrated by a time-split method with the thirdorder Runge-Kutta scheme (Skamarock et al., 2008; Kang and Cheong, 2017). The hyper-viscosity is applied to the predicted variables after one cycle of the Runge-Kutta time stepping is finished. When the hyper-viscosity is used as a viscosity, the built-in viscosity in the WRF model is turned off. The WRF model incorporates the second-order horizontal mixing and the sixth-order diffusion (that is, Laplacian operator and the third-order Laplacian operator) in an explicit manner. The sixth-order diffusion is applied as following:

where *Q* means any prognostic variable in the model except for the humidity and the surface pressure, and *α* is the diffusion coefficient expressed by *α* = 2^{−6}*p*^{−1} Δ*t*^{−1}*ε* with *ε* [=0.12] the filtering coefficient and *p*[=2] the number of passes of the diffusion scheme, respectively (Skamarock et al., 2008; Knievel et al., 2007). The filtering coefficient implies the amount of diffusion applied to the disturbance of twice the grid interval (2Δ*x*), and hence it should be given to lie between 0 and 1. In order to avoid Gibbs phenomenon due to the sixth-order diffusion, a flux limiter is introduced to the discretization in the WRF model (Knievel et al., 2007). While this prevents the monotonic variable having negative value, it is also known to reduce the effective diffusion at the small scales near the grid-size.

The tropical cyclone Chaba (18^{th} 2016) was simulated to demonstrate the performance of the highorder filter designed in this study. Chaba was generated at 0300 UTC 28 SEP 2016 and was reported to have turned into an extratropical cyclone at 0000 UTC 06 OCT 2016 by the Korean meteorological administration (KMA). In the mature stage, Chaba was classified as a medium-sized strong typhoon. The initial surface pressure of the tropical cyclone, obtained through a bogussing procedure, is presented in Fig. 5 in comparison with those without initialization. It should be noted that the tropical cyclone after initialization is represented with a more realistic structure than the tropical cyclone in the analysis data (Kurihara et al., 1993; Kwon and Cheong, 2010).

Time integration of the model was carried out for 3 days with or without the bogussing method. Fig. 6 shows the tracks of Chaba for four simulation cases together with the best track. All the simulations shown here presented the tracks quite close to the observation, and on average the track error becomes large with time as is expected. The difference among the simulations with different diffusion scheme produced very similar results to one another, implying insignificant sensitivity to the viscosity scheme. On the other hand, the track error by 72 h appears to be the largest for the case without the bogussing method and the WRF diffusion. The track error by 72 h in this case was 375 km, while the error of other cases are less than 280 km. Fig. 7 presents the time variation of the minimum pressure (*p*_{min}; central pressure) during three days for the same cases as in Fig. 6. Firstly it is noted that the results are quite sensitive to the viscosity scheme and the presence of the bogussing procedure (Kwon and Cheong, 2010; Cheong et al., 2011): None of them appears to have predicted the central pressure accurately, particularly around day 1 where the pressure has reached the lowest value. Without the bogussing procedure, the simulation even starts with considerably higher pressure than the observation, and keeps larger value until day 2. In the case of the WRF diffusion the pressure tends to increase monotonously with time, although it shows very weak decrease in the first 6 hours. When the diffusion is not used, however, the pressure drops with unrealistic large slope in the first 6 hours and shows increasing tendency as well as large deviation from the observation thereafter. In the case of the highorder filter, the simulation also looks to fail in reproducing reasonable pressure change around the period of 12-36 h. However, except this period, the predicted minimum pressure closely follows the observation. Comparing the results for the three days period, it may be stated that, on average, the simulation with the new filter provides better results than other cases.

One-day accumulated rainfall (05 OCT 2016) associated with the Chaba is presented in Fig. 8. The observation indicates that the rainband is distributed over the southern area of the Korean peninsula with heavy rainfall being found along the southern coastal line and over the Jeju Island. Four simulation results exhibited quite different rainfall distribution from the observation, which reflects the difficulty of accurate forecast of rainfall (e.g., Nguyen and Chen, 2014; Hogsett and Stewart, 2014). Among them, the result with the high-order filter seems to have reproduced better rainfall distribution than other three cases. The simulation with WRF diffusion shows too weak- and too much- rainfall over the southeastern and southwestern coast of the Korean peninsula, respectively. Furthermore, the rainfall was not simulated accurately over the central region of the Korean peninsula. It is quite interesting to point out that the simulation without the diffusion gave strongly biased distribution toward the South Sea. Without the bogussing procedure, the rainfall is forecast to spread over the much broader area than the observation. While strongly suggesting the necessity of the diffusion scheme in the model, the results clearly demonstrated the sensitivity to the details of the diffusion scheme.

## Summary and Conclusions

In this study a high-order filter using two-dimension finite element method with quadrilateral basis functions was developed to use in the numerical models on the spherical-surface limited area domain. The basis function consists of four shape-functions, each being defined on four grid-boxes around one shared gridpoint. The performance of the filter was evaluated in terms of the differential operators of the first order as well as the high-order Laplacian. It was shown that as a whole the accuracy is comparable to the filter based on the Fourier finite element method. The convergence rate of errors associated with the filter was also the same as the filter with Fourier finite element method. The high-order filter based on the two-dimension finite element method was incorporated in the WRF model as a hyper-viscosity. To compare the relative performances of the built-in hyper-viscosity and the new high-order filter, numerical simulations to predict the track and intensity of the tropical cyclone Chaba (18^{th} 2016) were carried out. It was found that the high-order filter performs better than the built-in viscosity scheme in that it provides the track and minimum pressure closer to the observation. It was also revealed that the error of the amount and spatial distribution of the rainfall was smaller for the highorder filter than the built-in viscosity scheme. The most important advantage of the filter developed in this study over the previous filter may be that the present filter does not need to use the Fourier transform because the calculation can be done using the gridpoint values rather than the spectral coefficients. This makes it much easier to implement the filter in numerical models while maintaining the same level of accuracy as those based on the spectral method in the zonal direction. Another advantage of the new filter may be that it can be implemented with ease and flexibility to the arbitrary grids. The new high-order filter deals with a huge matrix because the Fourier series is not incorporated in the zonal direction. The solution method is, therefore, very similar to that for cubed-sphere grid (Kang and Cheong, 2017) in many respect, one of the widely used non latitude-longitude grids. In a subsequent study, the parallel efficiency of the new filter will be discussed in depth, incorporating the details of Kang and Cheong (2017).