This page documents how to compute kinetic energy spectra from EAM atmosphere simulations.
Summary of best practices:
Spherical harmon transform available in NCL and pyNGL for lat/lon data.
Interpolate EAM native grid output to a (N+1)x2N lat/lon cap grid
Take N>=NE*6 (degrees of freedom, pole to pole)
Interpolate using TR’s “highorder” algorithm.
Need instantaneous output of (U,V) or (vor,div) on the GLL grid (not the PG2) grid.
For smoother results, need spectra from ~(how many?) snapshots
At NE256, interpolating and computing the spectra for each snapshot takes about 5min and 30GB of memory. At NE1024 this will start to become challenging.
KE Spectra
KE spectra is useful for understanding and tuning the dissipation mechanisms in the model. For a global model, it can be computed via a vector spherical harmonic transform of the instantaneous velocity field, or two spherical harmonic transforms of the instantaneous vorticity and divergence scalars. Spherical harmonics can be thought of as polynomials in Cartesian coordinates (x,y,z) restricted to the sphere. For each degree k there will be 2k+1 polynomials of total degree k, or 2k+1 spherical harmonics of degree n. To compute the KE spectra, E(k), we sum the coefficient squared over all degree n spherical harmonics.
The atmosphere is in geostrophic balance, resulting in a kinetic energy spectrum that scales like E(k) ~ k^-3 as k ranges from 3000km to 800km. But at 600km (wave number 134), the atmosphere transitions into a k^−5/3 scaling that holds down to a few km. The transition, or breakdown of the geostrophic balance, is still not fully explained. The −5/3 scaling suggests a transition to three-dimensional isotropic turbulence, but this explanation can be quickly ruled out. This was first observed by Nastrom and Gage, who suggested a combination of the known enstrophy generation from baroclinic instability with an unknown small scale energy source. Later investigators have speculated on where this source could come from: there are pure dynamical explanations and explanations which require convection and moist physics (such as thunderstorms).
Example KE spectra. Plotted is the “compensated” spectra, E(k)*k^5/3, to better illustrate the possible transition to ^-5/3. E(k) was computed from many flow snapshots and then E(k) is averaged over all the results to get the nice curves shown below. Unfortunately, one cannot compute E(k) from time averaged flow fields.
Computing KE Spectra
Interpolate to the correct lat/lon grid.
Computing the KE spectra requires performing spherical harmonic transforms. This can be done via the long standing SPHEREPACK Fortran package, which has an easy to use interface in NCL (and pyNGL?).
SPHEREPACK requires gridded data on either a Gaussian grid or an equal angle lat/lon “cap” grid which includes points at the poles. It relies on the FFT for much of the transforms, and for efficiency the grid should have dimensions N x 2N (Gaussian), or (N+1) x 2N (equal angle cap). The most common lat/lon data is an equal angle offset grid, which avoids pole points (see. CAM-FV Grid Overview ). Data on this offset grid would have to first be interpolated to a cap grid before computing spherical harmonic transforms.
Choice of interpolation algorithm and resolution
To compute spectral from native grid EAM output, we will first interpolate to a (N+1)x2N lat/lon cap grid. The KE spectra turns out to be sensitive to both the resolution of the lat/lon grid, and the algorithm used for interpolation. To determine sufficient resolution and the best algorithm, we use NE=256 (13km) Aqua planet simulation. This grid has 3072 points at the Equator, meaning it can potentially resolve spherical harmonics up to degree 1536 (corresponding to 2dx).
To look at the impacts of interpolation algorithm, we compute the spectra from a single snapshot using bilinear, integrated bilinear and TR’s highorder (integrated native SE shapefunctions), and two different resolutions. We plot E(k)/E_0(k), where we normalize all spectra by E_0, the spectra compute using the high order algorithm on an oversampled 3073x6144 lat/lon cap grid.
Three key points from the data below:
Computing E(k) at lower resolution or with certain algorithms introduces a roll off (reduced energy) and high frequencies. Since the behavior of the spectra, especially the rolloff starting around wave number 200 is of most interest, this should be avoided.
Using “PG2” output (cyan curve below) is not recommended - this is most likely due to the downscalings of the velocity when transforming from the dynamics to the PG2 physics grid, removing energy from small scales.
“intbilin” algorithm (green and brown) is not recommended (brow and green curves below). This algorithm is very nice for analysis, but does remove energy from small scales. “bilin” is a little better (orange, blue)
Best results are obtained by the highorder algorithm (red and straight line reference).
We thus recommend spectra be computed by first transforming to a lat/lon grid using TR’s high-order algorithm. We next look at the effects of interpolation grid resolution: Again plotting E(k)/E_0(k). Observations:
The HO algorithm is good enough (no signficant rolloff at wave number 200) at a resolution of 1537x3072.
The bilin algorithm requires almost twice as much resolution to get as good results.
Note that in the actual spectra plot, the differences are quite small, with the lower resolution spectra matching very well at low wave numbers, and only slightly below at high wave numbers.
Time Averaging
The spectra above from single snaphots have a lot of noise which can be removed by averaging. We next analyze how much averaging is necessary…