Open access peer-reviewed chapter

# Spectral Representation of Time and Physical Parameters in Numerical Weather Prediction

By Kristoffer Lindvall and Jan Scheffel

Submitted: April 16th 2018Reviewed: July 17th 2018Published: November 5th 2018

DOI: 10.5772/intechopen.80351

## Abstract

Numerical weather prediction (NWP) is a difficult task in chaotic dynamical regimes because of the strong sensitivity to initial conditions and physical parameters. As a result, high numerical accuracy is usually necessary. In this chapter, an accurate and efficient alternative to the traditional time stepping solution methods is presented; the time-spectral method. The generalized weighted residual method (GWRM) solves systems of nonlinear ODEs and PDEs using a spectral representation of time. Not being subject to CFL-like criteria, the GWRM typically employs time intervals two orders of magnitude larger than those of time-stepping methods. As an example, efficient solution of the chaotic Lorenz 1984 equations is demonstrated. The results indicate that the method has strong potential for NWP. Furthermore, employing spectral representations of physical parameters and initial values, families of solutions are obtained in a single computation. Thus, the GWRM is conveniently used for studies of system parameter dependency and initial condition error growth in NWP.

### Keywords

• NWP
• time-spectral
• chaotic
• error analysis
• initial condition

## 1. Introduction

Numerical methods are routinely used when modeling complex systems such as meteorological and atmospheric systems. These systems can be described by nonlinear ordinary and partial differential equations. Nonlinear systems in particular, as opposed to linear systems, constitute specific numerical challenges in terms of processing and memory requirements. In order to efficiently use the computer’s resources, it is important to choose a suitable numerical method [1].

Typically, when solving initial value ODEs, it is common practice to use finite difference methods (FDM) that discretize the temporal domain into finite steps. Information from the initial state and previously computed steps are used to estimate the solution at a “future” discrete time step. This is an intuitive and powerful method because we can easily visualize the cause and effect of each step. However, the problem with these techniques when solving complex systems is that small temporal steps are required in order for a solution to be found. In order to acquire a stable solution, the Courant-Friedrichs-Lewy (CFL) condition needs to be satisfied when employing more accurate, explicit FDM, thus limiting time step length.

Another approach is to use a time-spectral method, which takes a bird’s eye view of the problem [2, 3, 4, 5, 6, 7]. Instead of focusing on small local steps in time, an approximate solution that could fit the entire temporal domain is postulated. This approximate solution is set in the form of a finite series of basis functions. Thus, the unknowns of the equations are thereby changed from the physical variables to the coefficients of the solution ansatz. The procedure allows for highly accurate and global temporal solutions.

In the following, we will present application of the GWRM [7] to the Lorenz 1984 [8] differential equations. This is a set of three coupled, nonlinear, and chaotic ODEs, providing a basic model of the important process of Hadley tropical atmospheric circulation. In Section 2, Method, we present a short derivation of the GWRM. Section 3 will introduce the Lorenz equations and how they can be solved by the time-spectral method. In this section, we also show how the parameter dependence of the solution can be found from a single GWRM computation. We also show that, similarly, the dependence on initial conditions can be found from one computation only. A discussion is given in Section 4, followed by conclusion in Section 5.

## 2. Method

The time-spectral method employed here is the generalized weighted residual method (GWRM). The physical equations of this Galerkin method are projected onto a space of weighted orthogonal basis functions resulting in a set of algebraic equations. The residual of the differential equation is necessarily zero for the exact solution, thus the numerical method seeks to solve for the coefficients that minimize the residual in the interval where orthogonality holds.

In this procedure, all dimensions can be included, thus the term “generalized.” This means that the dependence on all temporal, spatial, and parameter domains can be included in the algebraic equations.

A key characteristic of Chebyshev polynomials, and a main reason for choosing them as basis functions, is their minimax property [9]. This property states that a Chebyshev series of order nis the best approximation of the same Chebyshev series of order n+1. The implication is that the highest mode can be neglected without deteriorating the approximate solution. This enables computation of nonlinear products entirely within spectral space. Actually, all GWRM numerical computations are carried out in spectral space. In contrast to, for example, Fourier series, Chebyshev series are also conveniently real. Any type of boundary condition can be employed, such as Dirichlet, Neumann, and periodic boundary conditions. Finally, Chebyshev polynomial series are known to converge rapidly when approximating smooth functions.

### 2.1. Generalized weighted residual method

Consider the following set of differential equations,

utFu=s,E1

where utxpis the unknown variable with a temporal, spatial, and parameter dependence. Fis a linear or nonlinear operator, and sis a source term. The GWRM solution ansatz takes the form, in the case of a singular physical dimension and one parameter,

utxpUτξP=k=0Kl=0Lm=0MaklmTkτTlξTmP.E2

Here, Tnx=cosnarccosxis the Chebyshev polynomial of the first kind, aklmare coefficients, and the prime (') denotes that the first term of each summation should be multiplied by 1/2. Since Chebyshev polynomials are orthogonal in the interval 11, variable transformations are applied for general intervals,

τ=tAtBt,ξ=xAxBx,P=pApBp,E3

where Az=z1+z0/2and Bz=z1z0/2, zbeing the variable t, x, and pwith their respective upper and lower bounds. The Picard integral is then applied to (1), after which the residual of the differential equation is formed,

R=utxput0xp+t0tFu+sdt.E4

The residual is multiplied by weighted Chebyshev polynomials and then integrated over the computational domain, whereafter it is set to zero. The equation takes the form

p0p1x0x1t0t1RTqτTrξTsPwtwxwpdtdxdp=0,E5

where the Chebyshev weight is wζ=1ζ21/2, ζbeing the transform variables τ, ξ, and P. Eq. (5) has the form of an integral minimum in calculus of variation, which states that if Ris smooth and differentiable, and the basis functions and weights are nonzero, then Rwill be minimized in the domain [10].

Continuing, we expand Eq. (5),

p0p1x0x1t0t1utxp[u(t0xp)+t0tFu+sdt]TqτTrξTsPwtwxwpdtdxdp=0,E6

yielding the following relations for each term of (6):

p0p1x0x1t0t1utxpTqτTrξTsPwtwxwpdtdxdp=BtBxBpπ23aqrsE7
p0p1x0x1t0t1ut0xpTqτTrξTsPwtwxwpdtdxdp=BtBxBpπ22πδq0brsE8
p0p1x0x1t0t1t0tFudt'TqτTrξTsPwtwxwpdtdxdp=BtBxBpπ23AqrsE9
p0p1x0x1t0t1t0tsdtTqτTrξTsPwtwxwpdtdxdp=BtBxBpπ23Sqrs.E10

where

t0tFudtk=0Kl=0Lm=0MAklmTkτTlξTmP,E11

and

t0tsdtk=0Kl=0Lm=0MSklmTkτTlξTmP.E12

Combining Eqs. (7a)(7d), we obtain the final set of algebraic equations

aqrs=2δq0brs+Aqrs+SqrsE13

where aqrsare the solution coefficients, δq0is the Kronecker delta function, brsrepresents the initial conditions, Aqrsis the spectral representation of the time-integrated linear or nonlinear operator F, and Sqrsrepresents the time-integrated source term. Eq. (10) is here defined for the truncated domain 0qK, whereas strictly the time integration renders K+1terms. For the spatial domain, we have 0rLBC, where LBC=LNBC, with NBCrepresenting the total number of boundary conditions for this spatial dimension. Boundary condition equations are thus added to (10) for problems including spatial derivatives. Here, we will be solving a set of time-dependent ODEs, thus only initial conditions are required to find a unique solution. Finally, it holds that 0sM.

Eq. (10) can be solved iteratively with a suitable root solver. Here, we use the semi-implicit root solver (SIR) [11] because it shows superior global stability compared to the Newton method (NM) with line-search while still maintaining a fast convergence.

For more details regarding GWRM basics, please see references [7, 12].

## 3. The GWRM, exemplified by the Lorenz 1984 equations

The Lorenz 1984 model is a set of three nonlinear ordinary differential equations that features chaotic behavior similar to that of meteorological systems [8]. It is consequently one of the simplest models of Hadley circulation due to the substantial amount of approximations postulated in order to arrive at this low order model. Thus, whereas the Lorenz 1984 model is not an accurate NWP model, it allows for a rigorous numerical analysis of simple, yet nonlinear chaotic behavior. This is the reason for employing them in the present work. The equations are

dXdt=Y2Z2αX+αQ,dYdt=XYβXZY+W,dZdt=βXY+XYZ.E14

For a more exhaustive description of the model, see [8]. Suffice it to say the variables X, Y, and Zrepresent certain meteorological systems such as wind currents and large-scale eddies. The coefficients α, β, Q, and Ware chosen within certain limits to act as damping, coupling, and amplification of the physical processes of the system.

Lorenz, in his 1984 article, posed the question, “What can such a simple model possibly tell us about the real atmosphere?”. His answer was that postulating certain hypothesis about the behavior of the real atmosphere, we can strengthen the reasoning behind the hypothesis qualitatively by the help of these models. However, that is not what we are concerned with here, but rather the model is used as a test to show how different numerical methods handle chaotic behavior. Furthermore, the model is useful to develop new techniques that can analyze such quantities as error growth, parameter uncertainties, and variable perturbations in a more efficient way.

### 3.1. Time intervals

The GWRM is a global method where a finite number of Chebyshev modes are used to represent the physics in a given temporal domain. Instead of using a large number of Chebyshev modes to obtain a single solution for the entire temporal domain, it is more efficient to split the domain into smaller subintervals, so that lesser Chebyshev modes can be used in each interval.

For controlling and maintaining the same numerical accuracy within each time interval, an automatic time-adaptive algorithm is employed. Time-adaptive techniques are easily implemented with Chebyshev polynomials. Since the absolute value of the Chebyshev coefficients decrease with increasing modes for well-resolved functions, the ratio of the absolute values of the highest modes to those of the lowest modes (see Eq. (12)) gives a direct estimate of how “good” the solution is. Thus, it holds that (for the example that (10) represents an ODE)

RK=aK+aK1a0+a10,KE15

where ak, k=0..K, are the Chebyshev coefficients with Kas the highest mode. By requiring that RK<ϵ, where εis a stipulated accuracy, the time subintervals to achieve this accuracy are automatically computed. We have found it efficient to use a balanced algorithm that tries to expand the time interval at about every 10th time interval by a factor 1.5, and that halves the interval whenever ϵaccuracy could not be achieved.

An interesting and useful property of Chebyshev polynomials is Tk1=1, so that

ut1k=0KakTk1=k=0Kak.E16

Thus, the exact solution at the end state ut1can be approximated by a sum of the Chebyshev coefficients. Eq. (13) conveniently allows for the previous end conditions, represented as Chebyshev coefficients, to be applied as initial conditions for the next interval. All computations of a series of connected intervals are thus calculated in spectral space. A GWRM time-interval pseudocode is provided in Appendix A. Only one domain, the temporal, is presented for the sake of clarity.

The pseudocode shows the following: the algebraic equations are collected in a vector ϕ, along with the initial conditions in the first element of the vector. The time domain is divided into subintervals in which the ϕequations are solved by SIR. The initial condition is equal to the end-state of the previous solution, see (13). If the convergence parameter “conv” of the solution is larger than the pre-set minimum error, then the current subinterval is halved and the ϕequations are solved by SIR. On the other hand, if the “conv” convergence criterion is satisfied, the algorithm immediately proceeds to solve the next time interval, the length of which is increased by a factor 1.5 typically every 10th interval. This procedure is done until the entire temporal domain is solved for. When all coefficients are solved for, the semianalytical solution ansatz is easily computed.

### 3.2. Parameter dependency

Instruments are continually measuring the weather’s temperature, wind speeds, and topological height and width to name a few parameters. However, the measurement devices cannot, at present, give a global coverage. Thus, meteorological parameters are approximate and interpolated.

The meteorological parameters are represented in the Lorenz model by α, β, Q, and W. The uncertainty of the measured parameters will inevitably lower the confidence of the solution, specifically important in weather prediction.

When studying parameter dependencies, the GWRM enables an interesting solution. Instead of assigning the parameter a single value, we can introduce a spectrum of values. Thus, we introduce a parameter dimension. This technique was first presented in [7], where it was applied to the viscosity parameter in the 1D nonlinear Burger equation. Later, it was applied to a two-dimensional magnetohydrodynamic problem in the paper by Riva et al. [13].

For the Lorenz ODEs, where we have temporal and parameter dependencies, the solution ansatz then takes the form,

utαUτP=k=0Km=0MakmTkτTmP.E17

In Eq. (14), a single parameter has been included; however, the procedure can be expanded to any number of parameters. The impact of the parameter αfrom the Lorenz equation (11) on the Xvariable is demonstrated in Figure 1. The result from a single GWRM computation is displayed, using the ansatz (14). It can be seen that varying αwith the other parameters of Eq. (11) held fixed, the solution is strongly dependent on αfor longer times. For lower values of α, the solution becomes more stable and fluctuates with a higher frequency than for higher values of α.

It is of interest to provide a quantitative measure of the effect of αon the solution. One approach would be to compare the deviations from a “base” run for different values of α. For simplicity, however, we have chosen to monitor the standard deviations as a measure of the deflection from the average value in the entire parameter interval. Standard deviations can be computed from the semianalytical solutions with the formula

σ=i=1NUiU¯2N1,E18

where Uiis the GWRM solution for a specific value of the parameter. The evaluations have been carried out at Nparameter points iat a time of interest. The parameter interval averaged solution is denoted by U¯. A table is presented below where comparisons are made with a base run with specified parameter values α, β, Q, and W. The solution at t=1.0is used. Average values U¯and standard deviations σ, representing the entire parameter interval are also given.

The standard deviations are supplied with indices to indicate the chosen parameter and variable for analysis. For example, σα,xstates that αis the parameter domain and the standard deviation of Xis being analyzed. One parameter at a time has been varied by ±20%, while the others are kept constant.

Parameter dependency example, using a single GWRM computation t=1.0, K=8, M=8
GWRM sol.X1P=1.311Y1P=0.046Z1P=1.541
α=0.2,0.3, α=0.25
U¯±σX¯±σα,xY¯±σα,yZ¯±σα,z
1.307±0.0790.061±0.3531.498±0.073
β=3.2,4.8, β=4.0
U¯±σX¯±σβ,xY¯±σβ,yZ¯±σβ,z
1.322±0.0680.036±0.5841.401±0.128
Q=6.4,9.6, Q=8.0
U¯±σX¯±σQ,xY¯±σQ,yZ¯±σQ,z
1.308±0.0960.071±0.4141.482±0.087
W=0.8,1.2, W=1.0
U¯±σX¯±σW,xY¯±σW,yZ¯±σW,z
1.310±0.0100.044±0.1101.536±0.015

It is immediately seen that the Yvariable is strongly affected by changes in all parameter values, whereas the Xand Zsolutions are more robust.

### 3.3. Initial condition uncertainty

Meteorological models generally include time-dependent PDEs and ODEs that require initial conditions as starting points. In mathematical terms, initial conditions are required to close the system of equations. If we are interested in computing a family of scenarios employing different initial conditions, the standard approach to determine the final states is to restart the computations from scratch, using a new initial condition each time.

Time-spectral methods offer a more convenient approach. Similarly, as when we computed parameter dependency in a single GWRM computation in the previous section, we will now study the effect of a spectrum of initial values on the solution of Eq. (11), employing a single GWRM computation.

Initial condition dependency example, using a single GWRM computation t=1.0, K=8, M=8
GWRM sol.X1P=1.311Y1P=0.046Z1P=1.541
X0P0.768,1.152, X=0.96
U¯±σX¯±σx,xY¯±σx,yZ¯±σx,z
1.306±0.0790.047±0.2051.521±0.087
Y0P1.320.88, Y=1.1
U¯±σX¯±σy,xY¯±σy,yZ¯±σy,z
1.319±0.1560.021±0.4401.460±0.092
Z0P0.4,0.6, Z=0.5
U¯±σX¯±σz,xY¯±σz,yZ¯±σz,z
1.309±0.0080.041±0.1321.535±0.008

An interesting, albeit challenging, feature of meteorological models is the inherent uncertainty of the initial conditions. As a result, how certain can we be of our numerical results? To be able to predict a chaotic system with absolute precision, one would need infinitely accurate measuring devices, computer, and numerical method. Since none of these are granted, error growth analysis is important to gauge how far into the future we can be confident of our predictions.

A classical error growth analysis employed for the Lorenz equations (11) would be the following. A base scenario Uk,1(where "k"denotes the variable) with initial conditions v=<X0,Y0,Z0>=<0.96,1.10,0.50>with a time window of t050is solved for. Subsequently, a number of perturbed scenarios Nare solved, denoted with superscript , where the initial conditions are perturbed an amount δ, taking values <0.001, for example. The error growth is then computed with the formula,

Ekt=1Nn=2NUk,nUk,12,k=1,2,and3.E19

Thus, in this traditional analysis, a large number of computations are needed, where the ODEs are solved for different perturbed initial conditions.

We suggest another approach. By use of an ansatz of the type (14), a spectrum of initial conditions is allowed in a single computation. In the table above, results are presented for three cases where Eq. (11) has been solved to time t=1.0. For the three cases, the base initial conditions X, Y, and Zas well as intervals for X0P, Y0P, and Z0Pare shown. Also, results for the base initial conditions at t=1are provided. For the solutions, where the initial conditions are allowed to vary in an interval, both averages over the interval are shown as well as corresponding standard deviations. The analysis shows that the Yvariable value at t=1is strongly dependent on the initial condition. In Figure 2, the time and initial condition dependence of the Xvariable is shown in a 3D diagram. Here, the run time has been extended to t=5.

It would also be of interest to study the effect of allowing the initial conditions for all variables X, Y, and Zto vary simultaneously. In order to accomplish this in a single GWRM run, the natural thing to do would be to extend the ansatz (14) to include three parameters, corresponding to the different initial conditions. Instead, we have chosen a simpler approach. We let X0δ=X+δ, Y0δ=Y+δ, and Z0δ=Z+δ, where the perturbation δis a single parameter that is applied to all three variables. An example computation is shown in Figure 3. The Xvariable here features a quite a different time evolution than in Figure 2, where only the Xvariable initial condition was perturbed.

In the table below, initial condition interval average values, as well as standard deviations, are provided for all the variables X, Y, and Zfor a run extending to t=20.

Initial condition dependency example, using a single GWRM computation t=20, K=8, M=8, δ<0.001
GWRM sol.X200=1.779Y200=0.483Z200=0.167
U¯±σX¯±σxY¯±σyZ¯±σz
1.787±0.0630.460±0.0620.099±0.226

Returning to error growth analysis as suggested by Eq. (16), reliable results can be computed if a sufficiently large domain has been spanned by the perturbed initial conditions. We have compared the classical approach of perturbed iterations (PI) with the new initial condition parameter dependency (ICPD) technique in Figure 4a and b. Of particular interest is the potential gain in CPU time, using the ICPD.

It is found that the ICPD technique with K=8and M=4achieves the same result of error growth as the PI scheme. The ICPD scheme computed the error growth in tCPU=2.7minas compared to the PI scheme which required tCPU=12.7min. Thus, a near fivefold increase in efficiency was obtained for this relatively simple case.

Numerically, Figure 4a is based on N=250runs with the PI scheme, whereas Figure 4b was computed in a single run, where N=1000different initial condition values was used.

## 4. Discussion

Is it more efficient to solve a system of differential equations using a spectral representation of the parameter domain instead of solving the system multiple times with different parameter values? Some points related to this discussion are that for the ICPD approach, using the GWRM,

• the parameter dependencies can easily be analyzed since GWRM solutions are analytical

• the time-spectral method is a high-order method, leading to high accuracy solutions

• all parameters in an interval are included, whereas certain regions of critical parameter dependence could be missed out with traditional PI methods

• the parameter domain can be split into intervals, that can be solved for separately, potentially spanning a larger parameter space more efficiently

In this chapter, a univariate analysis has been performed; however, a multivariate analysis can also be implemented in the same manor. It should be noted, however, that the more parameters introduced in the analysis, the more Chebyshev modes need to be solved for, which results in larger matrices and memory demands.

It is also of interest to address the accuracy of the present computations. The ICPD with K=8and M=4did not achieve high accuracy of the exact solution at the end time t=50as can be seen in Figure 5. To be more precise, the solution at higher times (t>25) does not represent the real dynamics because the introduced error has grown to an extent where the real solution is lost (see Figure 5). For long run times, the effect of uncertain initial conditions is muddled by the effect of numerical inaccuracy. On top of this, the introduction of the parameter dimension itself in the ansatz (14) is a source for inaccuracy, since the GWRM algorithm has to handle a larger set of Chebyshev polynomials in this case. The argument could be made that if a lower order (e.g., fourth order) solution ansatz cannot accurately represent the parameter “physics,” then the solution exhibits no predictive power if the prediction horizon has been exceeded. How could then the error growth be more accurately represented? The parameter Chebyshev coefficient could perhaps be increased, or the parameter domain could be split into intervals. These are interesting questions for future studies of ICPD methods related to NWP modeling.

## 5. Conclusion

Spectral methods have a long history when applied to the spatial domain in PDE modeling. The time-spectral method GWRM, demonstrated here, provides similar accuracy for the temporal domain. High-order methods, like the GWRM, may achieve resolutions much higher than that of finite difference methods for similar amount of work. Furthermore, when applied to meteorological systems, the time-spectral method can accurately and efficiently compute the physics in all spatial, temporal, and parameter domains.

A recent study, comparing the time-spectral method GWRM with commonly used time-stepping methods such as Runge-Kutta and other, high-order implicit methods, was carried out in [1]. Since then, the GWRM numerical algorithms have been further streamlined; see for example [14]. Moreover, it has been found that the Jacobian of the algebraic system of equations (10) need only be computed once by including the time interval length analytically. Thus, during the computations, the new time interval length can be substituted into the Jacobian before use. This decreases the computation time from tCPU=17.9minto tCPU=12.7min, roughly 30%, for 250perturbed scenarios with K=8and t050.

We have, in this chapter, also presented a new approach to error analysis using the initial condition parameter dependency (ICPD) technique, that is, by spectral representation of the parameter domain. In a single computation, the same error growth was reproduced as for the classical case (where a large number of runs were carried out for different initial conditions), with a computational time of tCPU=2.7min. This amounts to a near fivefold gain in CPU time efficiency.

Finally, it was shown that the ICPD technique is also successfully applicable for determining the effect of perturbed physical variables on the solution. In one GWRM computation only, it was found that the solutions to the Lorenz 1984 equations sensitively depend on the system parameter α.

Algorithm 1. GWRM time Intervals.
1: procedure
2:       ϕAk#Algebraic equations
3:       ϕ0A0+2.0IC#Apply initial conditionIC
4:       tacc0
5:       j1
6:       if tacc>Timethem
7:            conv1.
8:            while conv>rel.errordo
9:                  if j>1them
10:                        ϕ0A0+2.0xj10/2+xj11++xj1K.
11:                  xjSIRϕ.
12:                  convxK+xK1/x0+x1.
13:                  if conv<rel.errorthem
14:                        akxk#Save solution
15:                  else
16:                        ΔtΔt/2
17:            tacctacc+Δt
18:            if modpj10=0them
19:                   Δt1.5Δt
20:            jj+1

chapter PDF
Citations in RIS format
Citations in bibtex format

## More

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Kristoffer Lindvall and Jan Scheffel (November 5th 2018). Spectral Representation of Time and Physical Parameters in Numerical Weather Prediction, Understanding of Atmospheric Systems with Efficient Numerical Methods for Observation and Prediction, Lei-Ming Ma, Zhang Chang-Jiang and Feng Zhang, IntechOpen, DOI: 10.5772/intechopen.80351. Available from:

### Related Content

#### Understanding of Atmospheric Systems with Efficient Numerical Methods for Observation and Prediction

Edited by Lei-Ming Ma

Next chapter

By Feng Zhang, Yi-Ning Shi, Kun Wu, Jiangnan Li and Wenwen Li

First chapter

#### Introduction to Infrared Spectroscopy

By Theophile Theophanides

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

View all Books