27.1.3 An analysis of polar filtering

In the past, test case integrations for 5 years using a global domain with roughly $3^\circ$ resolution and realistic geometry, topography, initial conditions, and forcing yielded results which were different depending on whether fourier filtering or finite impulse filtering was used. In general, results using the finite impulse filter were inferior to those using the fourier filter based on a comparison with an integration using very small time steps and no filtering. For many years, the reason for this result was unknown and thought to be due to a code bug. However, the reason was not due to a code bug. The result can be explained with the following argument.

Assume time is discretized as $t=n\Delta\tau$ where n is the time step. At time t=0, assume that quantity T_n=0 is known. To predict (T_n=1), a change in T over the time step is computed as $\delta T_{0,1}$ (details are unimportant) and

$$T_{n=1} = T_{n=0} + \delta T_{0,1}$$ (27.23)

Then T_n=1 is filtered yielding

$$\tilde{T}_{n=1} = \tilde{T}_{n=0} + \tilde{\delta T_{0,1}}$$ (27.24)

To predict T_n=2, the next change in T is computed as $\delta T_{1,2}$ and

$$T_{n=2} = \tilde{T}_{n=1} + \delta T_{1,2}$$ (27.25)

After filtering, the result is

$$\tilde{T}_{n=2} = \tilde{\tilde{T}}_{n=1} + \tilde{\delta T_{1,2}}$$ (27.26)

which can be expanded to yield

$$\tilde{T}_{n=2} = \tilde{\tilde{T}}_{n=0} + \tilde{\tilde{\delta T_{0,1}}} + \tilde{\delta T_{1,2}}$$ (27.27)

Consider the term $\tilde{\tilde{T}}_{n=0}$. If a fourier filter is used with wavenumbers truncated above some critical cutoff, then two or more applications of this filter is the same as one application. The reason is that ``m'' applications of the filter is the same as multiplying ``m'' response functions together. However, since the response fuction of the simple finite impulse filter with weights (1/4, 1/2, 1/4) is a cosine, multiple applications yield successively more smoothing at lower wavenumbers.

An alternative to filtering prognostic variables is to filter the time tendencies. Using the above example, after $\delta T_{0,1}$ is computed, it can be filtered to yield $\tilde{\delta T_{0,1}}$. Then $\tilde{T}_{n=1}$ can be constructed as

$$\tilde{T}_{n=1} = T_{n=0} + \tilde{\delta T_{0,1}}$$ (27.28)

After computing and filtering to yield the next time tendency $\tilde{\delta T_{1,2}}$, the updated and filtered variable $\tilde{T}_{n=2}$ can be constructed as

$$\tilde{T}_{n=2} = \tilde{T}_{n=1} + \tilde{\delta T_{1,2}}$$ (27.29)

which can be expanded to yield

$$\tilde{T}_{n=2} = T_{n=0} + \tilde{\delta T_{0,1}} + \tilde{\delta T_{1,2}}$$ (27.30)

Comparing Equations 27.27 and 27.30, it can be seen that there are effectively fewer filter applications when time tendencies are filtered compared to when prognostic variables are filtered. Fourier filtering with a sharp cutoff gives essentially the same results using either method for the reason given above. However, the finite impulse filter will produce much less smoothing of long waves when time tendencies are filtered. Results using the explicit free surface and realistic forcing bear this out. For stability reasons within the explicit free surface scheme, it turns out that tracer tendencies and barotropic velocity tendencies can be filtered but baroclinic velocity must be filtered instead of the baroclinic tendency. Using the explicit free surface method and the above described filtering combination yields solutions which compare well with the unfiltered case. However, the filtering must be tunned to each geometric configuration.