MPAS | WRF | FDDA | Spectral nudging
Nudging
Nudging, also known as four-dimensional data assimilation (FDDA), is a method of keeping simulations close to analyses and/or observations or forcing data over the course of an integration by applying extra forcing terms to the model equations.
There are three types of nudging and some can be used in combination.
- Grid- or analysis-nudging simply forces the model simulation towards a series of analyses grid-point by grid-point.
- Observational- or station-nudging locally forces the simulation towards observational data.
- Spectral nudging only forces the model to the selected spectra of waves in the analysis.
These methods provide a four-dimensional analysis that is somewhat balanced dynamically, and in terms of continuity,while allowing for complex local topographical or convective variations. Such datasets can cover long periods, and have particular value in driving off-line air quality or atmospheric chemistry models.
WRF: Literature Review
The Weather Research and Forecasting (WRF) model's spectral nudging Four-Dimensional Data Assimilation (FDDA) system represents a sophisticated approach to constraining large-scale atmospheric circulation patterns while preserving mesoscale features within regional numerical weather prediction models. This technique, originally developed by von Storch et al. (2000) and later implemented in WRF by Miguez-Macho et al. (2004), offers a scale-selective method for incorporating observational constraints into dynamical downscaling applications.
Spectral nudging operates by decomposing the difference between model state variables and reference analysis fields into their spectral components through Fourier transforms, selectively nudging only the large-scale (long-wavelength) components while allowing smaller-scale features to develop freely within the model's physics. This approach represents a significant advancement over traditional analysis nudging techniques, which apply uniform corrections across all spatial scales without discrimination.
Spectral nudging operates by decomposing the difference between model state variables and reference analysis fields into their spectral components through Fourier transforms
- https://www2.mmm.ucar.edu/wrf/users/docs/technote/v4_technote.pdf
- Development and evaluation of spectral nudging strategy for the simulation of summer precipitation over the Tibetan Plateau using WRF (v4.0) | 2021
Literature Review
- Schubert-Frisius, Martina, et al. "Optimal spectral nudging for global dynamic downscaling." Monthly Weather Review 145.3 (2017): 909-927.
- Wedi, Nils P., Mats Hamrud, and George Mozdzynski. "A fast spherical harmonics transform for global NWP and climate models." Monthly Weather Review 141.10 (2013): 3450-3461.
- Park, Ja-Rin, Hyeong-Bin Cheong, and Hyun-Gyu Kang. "High-order spectral filter for the spherical-surface limited area." Monthly Weather Review 139.4 (2011): 1256-1278.
Spectral nudging
Spectral nudging is another way to nudge the model solution to either analysis or any forcing data. The key difference between this method and that of grid- or obs-nudging is that it is more selective in the scales one would like to nudge. In contrast to directly nudge model variables toward analysis or observation, this method decomposes the difference fields spectrally, and nudges toward the longer waves that correspond to the analysis. The spectral coefficients from the selected waves are then transformed back from the wave space to physical space and added to the model as tendencies. This technique was originally implemented to the Reginoal Atmospheric Modeling System, and later to ARW in 2010.
Spectral nudging uses the same input data as for grid-nudging, and nudges model variables u, v, temperature, and geopotential height. Water mixing ratio nudging was added in Version 4.0. Applications described above for grid-nudging can be considered using spectral nudging too.
Options to control the grid-nudging are also used for spectral nudging (such as nudging end time and ramping, nudging strength, controls of nudging in the boundary layer or lower levels and so on). In addition, there is an option to control the ramping in the vertical between the layer without nudging to the layer where the nudging is on. The number of waves to nudge in the x and y directions are user-defined.
module_fdda_spnudging.F
1 | |
spectral_nudging_filter_3dx
! Variables will stay in domain form since this routine is meaningless
! unless tile extent is the same as domain extent in E/W direction, i.e.,
! the processor has access to all grid points in E/W direction.
! There may be other ways of doing FFTs, but we haven't learned them yet...
1 | |
- fftpack5 routines
- call rfftmi(n,wsave,lensav,ier)
- forward transform: initialize coefficients, place in wsave
- call rfftmf( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
- do the forward transform
- filter all waves with wavenumber larger than nwave
- call rfftmb( lot, jump, n, inc, fin, lenr, wsave, lensav, work, lenwrk, ier )
- do the backward transform
- call rfftmi(n,wsave,lensav,ier)
FFTPACK5 provides a suite of routines for computing Fast Fourier Transforms (FFTs) of both real and complex data, including 1D and 2D transforms, as well as multiple vectors. These routines are designed for efficiency and flexibility, supporting various data types (real/complex, single/double precision) and lengths.
1 | |
WRF setting
- Spectral nudging with ERA5 | Jan 9, 2023
- In order to turn on spectral nudging, you will need to set
grid_fdda = 2. The wavenumber nudged is selected in namelist (e.g.,xwavenum,ywavenum, e.g. =3). Please choose the number so that- (domain size)/(wavenumber)=~1000 km in each direction.
- In order to turn on spectral nudging, you will need to set
- Gómez, B., and G. Miguez‐Macho. "The impact of wave number selection and spin‐up time in spectral nudging." Quarterly Journal of the Royal Meteorological Society 143.705 (2017): 1772-1786.
- when selecting smaller scales, the fine-scale contribution of the model is damped, thus making 1000km the appropriate scale threshold to nudge in order to balance both effects.
- \(119 \times 105\) horizontal points at 36km resolution and 33 vertical, \(119 \times 36 = 4284km\)
- Silva, Natália Pillar da, and Ricardo de Camargo. "Impact of wave number choice in spectral nudging applications during a South Atlantic convergence zone event." Frontiers in Earth Science 6 (2018): 232.
- 1000km - R = the Rossby radius length for most mesoscale studies
ERA5 (Not the point, Ignored)
Spherical harmonic transforms
Performing a global spherical Fourier transform on ERA5 data in Python typically involves using spherical harmonic transforms, which are the generalization of the Fourier transform to the sphere. This is particularly relevant in the context of global weather forecasting and climate modeling, where the Earth's spherical geometry needs to be accurately represented.
- Q: how to define high wavenumber to filter Rossby waves?
pyshtools
SHTOOLS/pyshtools is a Fortran-95/Python library that can be used for spherical harmonic transforms, multitaper spectral analyses, expansions of gridded data into Slepian basis functions, standard operations on global gravitational and magnetic field data.
1 | |
- Tutorials & guides
- The easiest way to start learning about pyshtools is to see it in action in a Jupyter notebook.
- https://shtools.github.io/SHTOOLS/python-examples.html
Spherical Harmonic Transform Libraries:
- Utilize Python libraries designed for spherical harmonic transforms. pyshtools is a prominent library for this purpose, offering functions for forward and inverse spherical harmonic transforms.
- These libraries handle the complexities of representing functions on a sphere and performing the appropriate transforms.
1 | |
Empirical orthogonal function expansion (EOF) (Not used)
- have to find out needed components as Rossby waves (2D).
XY Fourier transform
real-space to k-space
\[\begin{align}f(x) &= \sum_{n=1}^{N} F(k_n) e^{ik_n x} \qquad x \in [0,L] \\\text{Assume} \qquad f(x \pm L) &= f(x) \\\Rightarrow e^{\pm i k_n L} &= 1 \qquad \forall n \\& k_n (n = 0,1,...,N)\\k_n L &= 2n \pi \qquad m \in \mathbb{Z} \\k_n &= \frac{2 \pi}{L} n \\dk &= \frac{2 \pi}{L} \\\because \lambda &= \frac{2 \pi}{k} \\\lambda_n &= \frac{L}{n}\\\end{align}\tag{1}\]
| \(n\) | \(k_n\) | \(\lambda\) |
|---|---|---|
| 0 | \(k_0 = 0\) | \(\infty\) |
| 1 | \(k_1= \frac{2\pi}{L}\) | \(L\) |
| 2 | \(k_1= \frac{2\pi}{L} \times 2\) | \(\frac{L}{2}\) |
| 3 | \(k_1= \frac{2\pi}{L} \times 3\) | \(\frac{L}{3}\) |
| \(\cdots\) | \(\cdots\) | \(\cdots\) |
| N | \(k_1= \frac{2\pi}{L} \times N\) | \(\frac{L}{N}\) |
np.fft2
Performing 2D Fourier transforms (forward and backward) in Python is typically done using the numpy.fft module, particularly the fft2 and ifft2 functions
<div class="fold"> <div class="fold-title fold-info collapsed" data-toggle="collapse" href="#collapse-90469042" role="button" aria-expanded="false" aria-controls="collapse-90469042"> <div class="fold-arrow">▶</div>FFT </div> <div class="fold-collapse collapse" id="collapse-90469042"> <div class="fold-content"> 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
import numpy as np
import matplotlib.pyplot as plt
# Create a sample 2D image (e.g., a 2D Gaussian)
# L = [0, 1]
def create_sample_image(size=1280):
x = np.linspace(0.0, 1.0, size)
y = np.linspace(0.0, 1.0, size)
X, Y = np.meshgrid(x, y)
image = np.sin(2*X*2*np.pi) + 0.5*np.sin(3*X*2*np.pi) + 0.2*np.sin(10*X*2*np.pi) + 0.1*np.sin(25*X*2*np.pi)
return image
# Perform 2D Fourier Transform, filter high wavenumbers, and inverse transform
def fourier_filter(image, Nx_kmax, Ny_kmax):
# Step 1: Perform 2D FFT
f_transform = np.fft.fft2(image)
f_transform_shifted = np.fft.fftshift(f_transform) # Shift zero to center
# Step 2: Create a low-pass filter mask
# Create wavenumber grid
N = image.shape[0] # Grid size
print("N = ", N)
# L = 20.0 # Domain length
# dx = L / N
# dk = 2 * np.pi / L
# k = np.linspace(-N*dk/2, N*dk/2, N, endpoint=False)
k = np.linspace(-N/2, N/2, N, endpoint=False)
KX, KY = np.meshgrid(k, k)
# Define high-wavenumber filter (e.g., remove wavenumbers above N_kmax)
filter_mask_KX = np.abs(KX) <= Nx_kmax
filter_mask_KY = np.abs(KY) <= Ny_kmax
# Step 3: Apply the filter in k domain
f_transform_filtered = f_transform_shifted * filter_mask_KX
f_transform_filtered = f_transform_filtered * filter_mask_KY
# Step 4: Inverse FFT
f_transform_unshifted = np.fft.ifftshift(f_transform_filtered) # Shift back
filtered_image = np.fft.ifft2(f_transform_unshifted)
filtered_image = np.real(filtered_image) # Take magnitude to ensure real output
return filtered_image, f_transform, f_transform_filtered
# Main execution
if __name__ == "__main__":
# Generate sample image
size = 1280
x = np.linspace(0.0, 1.0, size)
k = np.linspace(-size/2, size/2, size, endpoint=False) # unit = (2*np.pi/1.0)
# m_lambda = np.linspace(-size/2, size/2, size, endpoint=False) # unit = L/
image = create_sample_image(size)
# Apply Fourier transform and filtering
Nx_kmax = 6 # size/2
Ny_kmax = 30
filtered_image, f_transform, f_transform_filtered = fourier_filter(image, Nx_kmax, Ny_kmax)
#===============================================================
# One Dim
plt.figure(figsize=(12, 4))
plt.suptitle(f"Cut-off wavenumber \n Nx_kmax = {Nx_kmax}", fontsize=24)
plt.subplot(131)
plt.imshow(image, cmap='gray')
plt.title('Original Image')
plt.colorbar()
plt.subplot(132)
plt.imshow(np.log(np.abs(f_transform) + 1), cmap='gray')
plt.title('Fourier Transform (log scale)')
plt.colorbar()
plt.subplot(133)
plt.imshow(filtered_image, cmap='gray')
plt.title('Filtered Image')
plt.colorbar()
plt.tight_layout()
plt.show()
plt.close()
#===============================================================
# Visualize results
plt.figure(figsize=(12, 8))
plt.subplot(221)
plt.plot(x, image[0, :])
plt.grid()
plt.title('Original Image')
plt.subplot(222)
plt.plot(k, np.log(np.abs(f_transform[0, :]) + 1), 'ro')
plt.xlabel('k [unit = (2*np.pi/L)]')
plt.grid()
plt.title('Fourier Transform (log scale)')
plt.subplot(223)
plt.plot(x, filtered_image[0, :], label='Filtered Image')
plt.title('Filtered Image')
plt.grid()
plt.legend()
plt.subplot(224)
plt.plot(x, image[0, :], label='Original Image')
plt.plot(x, filtered_image[0, :], label='Filtered Image')
plt.title('Filtered Image')
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
plt.close()
#===============================================================
# k and \lambda : Fourier Transform (log scale)
plt.figure(figsize=(12, 4))
plt.subplot(121)
plt.plot(k, np.log(np.abs(f_transform[0, :]) + 1), 'ro')
plt.xlabel('k [unit = (2*np.pi/L)]')
plt.grid()
plt.title('Fourier Transform (log scale)')
plt.subplot(122)
plt.plot(k[:10], np.log(np.abs(f_transform[0, :10]) + 1), 'ro')
# plt.xlabel('m [unit = (L/\lambda)]')
plt.xlabel('k [unit = (2*np.pi/L)]')
plt.grid()
plt.title('Fourier Transform (log scale)')
plt.tight_layout()
plt.show()
plt.close()
#===============================================================
# k and \lambda : Fourier Transform
plt.figure(figsize=(12, 4))
plt.subplot(121)
plt.plot(k+size/2, np.abs(f_transform[0, :]), 'ro')
plt.xlabel('k [unit = (2*np.pi/L)]')
plt.grid()
plt.title('Fourier Transform')
plt.subplot(122)
plt.plot(k[:30]+size/2, np.abs(f_transform[0, :30]), 'ro')
# plt.xlabel('m [unit = (L/\lambda)]')
plt.xlabel('k [unit = (2*np.pi/L)]')
plt.grid()
plt.title('Fourier Transform')
plt.tight_layout()
plt.show()
plt.close()
#===============================================================
plt.plot(filtered_image[0, :]- image[0, :], 'ro')
plt.title('Filtered Image - Original Image')
plt.grid()
plt.tight_layout()
plt.show()
</div> </div></div>Checking
Orignal signal
Add some noises
More wavenumbers
2D waves
1 | |
- plot at
y = 0.25
Results
Thinking
- Rationale for the Approach
- We can express the relationship as (_m = ). If the length (L) decreases while maintaining a constant mass (m), the resulting wavelength (_m) becomes shorter, indicating that the corresponding wave number (k) increases.
This increase in (k) shifts the distribution towards higher values, resulting in a greater number of higher (k) values being included in the analysis.
Apply on MPAS-A/ERA5
- The calculation of spectral nudging is inside the source code.
- The MPAS-A is unstructed mesh, how to do FFT?
- Global mesh --> cartesian or spherical harmonic function
- unstructed mesh --> regrid to regular mesh?
- if yes, delta = model - ERA5 --> regrid --> FFT --> cutoff --> inverse FFT --> regrid back --> add filtered_delta ?