NNMD

class fridom.framework.projection.nnmd.NNMD(mset: ModelSettingsBase, order=3, epsilon=None, use_model=True, time_step_factor=0.01, use_discrete=True, enable_dealiasing=True)

Bases: Projection

Nonlinear normal mode decomposition

Parameters

msetModelSettings

The model settings.

orderint

The order of the balanced state.

epsilonfloat (default: None)

The small parameter epsilon.

use_discretebool (default: True)

Whether to use the discrete eigenvectors vectors.

enable_dealiasingbool (default: True)

Whether to enable dealiasing when computing the nonlinear term.

Description

Nonlinear normal mode decomposition (NNMD) is a method to filter out all fast waves from a given state \(\boldsymbol z\). The obtained state \(\boldsymbol z_b\) is called a balanced state. For a linear system, such a balancing operation would be the projection onto the eigenspace with the smallest eigenvalue, e.g. the geostrophic mode for ocean models. In contrast to the simple projection onto the linear geostrophic mode, NNMD takes the nonlinear terms into account.

The key concept of NNMD emerged independently through the works of Machenhauer (1977) [1] and Bear & Tribbia (1977) [2] . Building on this foundation, Warn et al. (1995) [3] further developed this approach by expanding the fast linear normal modes in a power series expansion in terms of a small parameter, the Rossby number. Here we use a small modification of the Warn et al. (1995) [3] method, that is described in more detail by Eden et al. (2019) [4].

Requirements

Form of the System

To apply nonlinear normal mode decomposition, we require a system of equations that can be written in spectral space as:

\[\partial_t \boldsymbol z = -i \mathbf A \cdot \boldsymbol z + \epsilon \boldsymbol N (\boldsymbol z)\]

where \(\boldsymbol z(\boldsymbol k, t)\) is the state vector, \(\boldsymbol k\) is the wave vector, \(t\) is the time, \(\mathbf A\) is a matrix, \(\epsilon\) is a small parameter and \(\boldsymbol N(\boldsymbol z)\) is a nonlinear term. The \(j\)-th. component of the nonlinear term is given by:

\[N_j(\boldsymbol z) = \boldsymbol z * (\mathbf G_j \cdot \boldsymbol z)\]

where the star \(*\) denotes a convolution, and \(\mathbf G_j\) is a matrix. The nonlinear term \(\boldsymbol N\) typically corresponds to the advection term, and \(\mathbf G_j \cdot \boldsymbol z\) to the gradient of the \(j\)-th component of \(\boldsymbol z\).

Scaling of the system

The system must be scaled, e.g. the magnitude of \(\mathbf A\), and of the nonlinear term \(\boldsymbol N\) should be of order one.

Eigenvectors and Eigenvalues

The eigenvectors \(\boldsymbol q_j\) and eigenvalues \(\lambda_j\) of the linear system matrix \(\mathbf A\) should be known. Further, projection vectors \(\boldsymbol p_j:math:\) are required, such that

\[\boldsymbol p_i \cdot \boldsymbol q_j = \delta_{i,j}\]

with the Kronecker-Delta symbol \(\delta_{i,j}\). The first eigenvalue \(\lambda_0\) must be zero.

Calculation of the balanced state

In this section we present the equations from which the balance state is computed. For more details on how these equations are derived and why the obtained state is balanced see the derivation section.

Nonlinear normal mode decomposition is based around the idea on expanding a state around a small parameter \(\epsilon\) and balancing the system of equations in each order of \(\epsilon\). The full balanced state \(\boldsymbol z_b\) up to order \(N\) is then given by:

\[\boldsymbol z_b = z_{0,0,0}~\boldsymbol q_0 + \sum_{n=1}^N \epsilon^n \sum_{j={1,2}} z_{j,n,0}~\boldsymbol q_j \]

where \(N\) is the maximum order to which the state is balanced, and \(z_{j,n,k}\) are amplitudes, where - \(j\) correspond to the \(j\)-th linear normal mode - \(n\) correspond to the \(n\)-th order of the power series expansion - \(k\) correspond to the \(k\)-th partial time derivative For the slow linear normal modes (\(j=0\)), these amplitudes are given by:

\[\begin{split}z_{0,0,k} = \begin{cases} \boldsymbol p_0 \cdot \boldsymbol z & \text{for} & k=0\\ \boldsymbol p_0 \cdot \boldsymbol I_0^{(k-1)} &\text{else} \end{cases}\end{split}\]

where \(\boldsymbol z\) is the state to be balanced, and \(\boldsymbol I_n^{(k)}\) are interaction terms that are given below. For the two fast linear normal modes \(j=1,2\), the amplitudes are given by:

\[\begin{split}z_{j,n,k} = \begin{cases} 0 & \text{for} & n=0\\ \frac{i}{\lambda_j} \left( z_{j,n-1,k+1} - \boldsymbol p_j \cdot \boldsymbol I_{n-1}^{(k)}\right) &\text{else} \end{cases}\end{split}\]

For \(n=0\), the interaction terms \(\boldsymbol I_n^{(k)}\) are given by:

\[\boldsymbol I_0^{(k)} = \sum_{m=0}^{k} \binom{k}{m} \boldsymbol S \left( z_{0,0,k-m} ~ \boldsymbol q_0 , z_{0,0,m} ~ \boldsymbol q_0 \right)\]

and for \(n>0\):

\[\begin{split}\begin{eqnarray} \boldsymbol I_n^k = \sum_{m=0}^k \begin{pmatrix}k \\ m \end{pmatrix} \left [ 2 \boldsymbol S \left( z_{0,0,k-m} ~ \boldsymbol q_0, \sum_{j=1,2} z_{j,n,m} ~ \boldsymbol q_j \right) \right. \\ \left. + \sum_{i=1}^{n} \boldsymbol S \left( \sum_{j=1,2} z_{j,i,k-m} ~ \boldsymbol q_j , \sum_{j=1,2} z_{j,n-i,m} ~ \boldsymbol q_j \right) \right] \end{eqnarray}\end{split}\]

where \(\binom{k}{m}\) are binomial coefficients and \(\boldsymbol S\) is a symmetrical bilinear form, defined as:

\[\boldsymbol S(\boldsymbol z_1 , \boldsymbol z_2) = \frac{1}{2} \left [ \boldsymbol N (\boldsymbol z_1 + \boldsymbol z_2) - \boldsymbol N (\boldsymbol z_1) - \boldsymbol N(\boldsymbol z_2) \right]\]

Note that the amplitudes \(z_{j,n,k}\) can be straightforward calculated from these equations, since \(z_{j,n,k}\) only depends on terms of order \(n-1\) and smaller, and \(z_{j,0,k}\) can be calculated from terms with time derivative order \(k'\) smaller than \(k\).

Derivation

We start with the spectral system of equations and multiply it with the projection vector \(\boldsymbol p_j\) from the left:

\[\partial_t z_0 = \epsilon ~\boldsymbol p_0 \cdot \boldsymbol N(\boldsymbol z) \quad , \quad \partial_t z_j = -i \lambda_j z_j + \epsilon ~\boldsymbol p_j \cdot \boldsymbol N(\boldsymbol z)\]

where \(z_j\) is amplitude of the \(j\)-th normal mode:

\[z_j \equiv \boldsymbol p_j \cdot \boldsymbol z\]

Note that we consider the zero mode \(j=0\) separately, since the eigenvalue is zero. For a balanced state \(\boldsymbol z_b\), the tendency of the state should be small:

\[\frac{\mathcal O(\partial_t \boldsymbol z_b)}{\mathcal O(\boldsymbol z_b)} = \epsilon \]

Introducing a slow time scale \(T=\epsilon t\), the balance condition becomes:

\[\partial_t = \epsilon \partial_T \quad \implies \quad \frac{\mathcal O(\partial_T \boldsymbol z)}{\mathcal O(\boldsymbol z)} = 1\]

which is easier to handle. In the slow time scale, the system of equations becomes:

\[\partial_T z_0 = \boldsymbol N (\boldsymbol z) \quad , \quad \epsilon \partial_T z_j = -i \lambda_j z_j + \epsilon \boldsymbol N(\boldsymbol z)\]

these equations will hereafter be referred to as the slow time scale tendency equations. While these equations still represent the full system of equations and thus, the solution space is not yet limited to balanced states. However, together with the balance condition above, the solution space is restricted to slowly evolving states, i.e. balanced states. However to find the balanced state \(\boldsymbol z_b\) of a given state \(\boldsymbol z\), we need a boundary condition for solving these equations. Such a boundary condition could for example be \(\boldsymbol z - \boldsymbol z_b\) is minimal. However, this boundary condition would be hard to work with. A simpler boundary condition is to take the linear slow mode as a base point coordinate, i.e.:

\[\boldsymbol p_0 \cdot \boldsymbol z_b = \boldsymbol p_0 \cdot \boldsymbol z = z_0\]

With this boundary condition, \(z_0\) is given, and we can solve the above balance equations for \(z_j\). To do so, we express the fast modes as a power series expansion:

\[z_j = \sum_{n=0}^\infty \epsilon^n z_{j,n}\]

Such that the full balanced state is given by:

\[\boldsymbol z_b = z_{0}~\boldsymbol q_0 + \sum_{n=0}^\infty \epsilon^n \sum_{j={1,2}} z_{j,n}~\boldsymbol q_j \]

We now insert this power series expansion into the slow time scale tendency equations, and sort the terms to the order of \(\epsilon\). We then combine the obtained equations with the balance condition by balancing the tendency equation in every order. We will get separate equations for the terms of order \(\epsilon^0\), of order \(\epsilon^1\), etc. By definition, the balance condition is satisfied when the terms balance in each order. However, before we insert the power series expansion, we need to find a way to deal with the nonlinear term.

Expansion of the Nonlinear Term

Using that the \(j\)-th component of the nonlinear term is given by \(\boldsymbol z * (\mathbf G_j \cdot \boldsymbol z)\) , we find:

\[\begin{split}\begin{eqnarray} N_j(\boldsymbol z_1 + \boldsymbol z_2) &=& \boldsymbol z_1 * (\mathbf G_j \cdot \boldsymbol z_1) + \boldsymbol z_1 * (\mathbf G_j \cdot \boldsymbol z_2) + \boldsymbol z_2 * (\mathbf G_j \cdot \boldsymbol z_1) + \boldsymbol z_2 * (\mathbf G_j \cdot \boldsymbol z_2) \\ &=& N_j(\boldsymbol z_1) + N_j(\boldsymbol z_2) + \boldsymbol z_1 * (\mathbf G_j \cdot \boldsymbol z_2) + \boldsymbol z_2 * (\mathbf G_j \cdot \boldsymbol z_1) \\ &\equiv& N_j(\boldsymbol z_1) + N_j(\boldsymbol z_2) + 2 S_j (z_1, z_2) \end{eqnarray}\end{split}\]

with the symmetrical bilinear form \(S_j\):

\[\begin{split}\begin{eqnarray} S_j (z_1, z_2) &=& \frac{1}{2} \left [ \boldsymbol z_1 * (\mathbf G_j \cdot \boldsymbol z_2) + \boldsymbol z_2 * (\mathbf G_j \cdot \boldsymbol z_1) \right] \\ &=& \frac{1}{2} \left [ N_j(\boldsymbol z_1 + \boldsymbol z_2) - N_j(\boldsymbol z_1) - N_j(\boldsymbol z_2) \right] \end{eqnarray}\end{split}\]

The first line shows that \(S_j\) is symmetric and linear in both arguments. The second line shows how to compute \(S_j\) from \(N_j\). Note that the nonlinear term \(\boldsymbol N\) can be expressed with the symmetrical bilinear form:

\[\boldsymbol N(\boldsymbol z) = \boldsymbol S(\boldsymbol z, \boldsymbol z)\]

Inserting the full balanced state and it’s power series expansion yields:

\[\begin{split}\begin{eqnarray} \boldsymbol S(\boldsymbol z_b, \boldsymbol z_b) &=& \boldsymbol S \left( \boldsymbol z_0 + \sum_{n=0}^\infty \epsilon^n \boldsymbol z_{f,n} ~,~ \boldsymbol z_0 + \sum_{i=0}^\infty \epsilon^i \boldsymbol z_{f,i} \right) \\ &=& \mathcal S(\boldsymbol z_0, \boldsymbol z_0) + \sum_i \epsilon^i \boldsymbol S(\boldsymbol z_0, \boldsymbol z_{f,i}) + \sum_n \epsilon^n \boldsymbol S(\boldsymbol z_{f,n}, \boldsymbol z_0) + \sum_{n,i} \epsilon^{n+i} \boldsymbol S(\boldsymbol z_{f,n}, \boldsymbol z_{f,i}) \\ &=& \mathcal S(\boldsymbol z_0, \boldsymbol z_0) + \sum_{n=0}^\infty \epsilon^n \left( 2 \boldsymbol S(\boldsymbol z_0 , \boldsymbol z_{f,n}) + \sum_{i=0}^n \boldsymbol S(\boldsymbol z_{f,i}, \boldsymbol z_{f,n-i}) \right) \\ &\equiv& \sum_{n=0}^\infty \epsilon^n \boldsymbol I_n \end{eqnarray}\end{split}\]

with

\[\boldsymbol z_0 \equiv z_0 ~ \boldsymbol q_0 \quad , \quad \boldsymbol z_{f,n} \equiv \sum_{j=1,2} z_{j,n} ~ \boldsymbol q_j\]

and the interaction term \(\boldsymbol I_n\):

\[\begin{split}\begin{eqnarray} \boldsymbol I_0 &=& \boldsymbol S(\boldsymbol z_0 + \boldsymbol z_{f,0} ~,~ \boldsymbol z_0 + \boldsymbol z_{f,0}) \\ \boldsymbol I_n &=& 2 \boldsymbol S(\boldsymbol z_0 , \boldsymbol z_{f,n}) + \sum_{i=0}^n \boldsymbol S(\boldsymbol z_{f,i}, \boldsymbol z_{f,n-i}) \quad \quad \text{for} \quad n>0 \end{eqnarray}\end{split}\]

Balance in Every Order

Inserting the power series in the slow time scale tendency equation of the fast wave modes \(j=1,2\) yields:

\[\begin{split}\begin{eqnarray} 0 &=& \epsilon ~ \partial_T z_j + i\lambda_j z_j - \epsilon ~ \boldsymbol p_j \cdot \boldsymbol S(\boldsymbol z_b, \boldsymbol z_b) \\ &=& \sum_{n=0}^\infty \epsilon^{n+1} \partial_T z_{j,n} + i \epsilon^n \lambda_j z_{j,n} - \epsilon^{n+1} \boldsymbol p_j \cdot \boldsymbol I_n \\ &=& i\lambda_j z_{j,0} + \sum_{n=1}^\infty \epsilon^{n} \left( \partial_T z_{j,n-1} + i \lambda_j z_{j,n} - \boldsymbol p_j \cdot \boldsymbol I_{n-1} \right) \end{eqnarray}\end{split}\]

balancing this in every order yields

\[\begin{split}\begin{eqnarray} z_{j,0} &=& 0 \\ z_{j,n} &=& \frac{i}{\lambda_j} \left( \partial_T z_{j,n-1} - \boldsymbol p_j \cdot \boldsymbol I_{n-1} \right) \quad \quad \text{for} \quad n>0 \end{eqnarray}\end{split}\]

The slow time derivative

For the calculation of the term \(z_{j,n}\) we require the slow time derivative of the term \(z_{j,n-1}\). We obtain an analytical expression for this term by taking the \(k\)-th slow time derivative of \(z_{j,n}\):

\[\begin{split}\begin{eqnarray} \partial_T^k z_{j,n} &=& \frac{i}{\lambda_j}(\partial_T^{k+1} z_{j,n-1} - \boldsymbol p_j \cdot \partial_T^k \boldsymbol I_{n-1}) \\ \Leftrightarrow z_{j,n,k} &=& \frac{i}{\lambda_j}(z_{j,n-1,k+1} - \boldsymbol p_j \cdot \boldsymbol I_{n-1}^{(k)}) \end{eqnarray}\end{split}\]

where the third index denotes the order of the derivative. To calculate the \(k\)-th derivative of the interaction term \(\boldsymbol I_n\), we need to take the derivative of the symmetrical bilinear form:

\[\begin{split}\begin{eqnarray} \partial_T S_j (\boldsymbol z',\boldsymbol z'') &=& \frac{1}{2}\partial_T \left( \boldsymbol z' * (\mathbf G_j \cdot \boldsymbol z'') + \boldsymbol z'' * (\mathbf G_j \cdot \boldsymbol z') \right) \\ &=& \frac{1}{2}\left( (\partial_T \boldsymbol z') * (\mathbf G_j \cdot \boldsymbol z'') + \boldsymbol z' * (\mathbf G_j \cdot \partial_T \boldsymbol z'') + (\partial_T \boldsymbol z'') * (\mathbf G_j \cdot \boldsymbol z') + \boldsymbol z'' * (\mathbf G_j \cdot \partial_T \boldsymbol z') \right) \\ &=& S_j(\partial_T \boldsymbol z', \boldsymbol z'') + S_j(\boldsymbol z', \partial_T \boldsymbol z'') \end{eqnarray}\end{split}\]

We search for a formula for the derivative of order \(k\). For this we take a look at the second order derivative:

\[\begin{split}\begin{eqnarray} \partial_T^2 \boldsymbol S(\boldsymbol z', \boldsymbol z'') &=& \partial_T \left[ \boldsymbol S(\partial_T \boldsymbol z', \boldsymbol z'') + \boldsymbol S(\boldsymbol z', \partial_T \boldsymbol z'') \right] \\ &=& \boldsymbol S(\partial_T^2 \boldsymbol z', \boldsymbol z'') + 2\boldsymbol S(\partial_T \boldsymbol z' , \partial_T \boldsymbol z'') + \boldsymbol S(\boldsymbol z', \partial_T^2 \boldsymbol z'') \end{eqnarray}\end{split}\]

This is the third row of pascal’s triangle. For the \(k\)-th. derivative, we find:

\[\partial_T^k \boldsymbol S(\boldsymbol z', \boldsymbol z'') = \sum_{m=0}^k \binom{k}{m} \boldsymbol S(\partial_T^{k-m} \boldsymbol z', \partial_T^m \boldsymbol z'')\]

where \(\binom{k}{m}\) is the binomial coefficient. Using this identity, the \(k\)-th derivative of the interaction term is given by:

\[\begin{split}\begin{eqnarray} \boldsymbol I_0^{(k)} &=& \sum_{m=0}^k \binom{k}{m} \boldsymbol S(\boldsymbol z_{0,0,k-m} , \boldsymbol z_{0,0,m}) \\ \boldsymbol I_n^{(k)} &=& \sum_{m=0}^k \binom{k}{m} \left[ 2 \boldsymbol S(\boldsymbol z_{0,0,k-m} , \boldsymbol z_{f,n,m}) + \sum_{i=0}^n \boldsymbol S(\boldsymbol z_{f,i,k-m}, \boldsymbol z_{f,n-i,m}) \right] \quad \quad \text{for} \quad n>0 \end{eqnarray}\end{split}\]

with \(\boldsymbol z_{0,0,k} = \partial_T^k \boldsymbol z_0\), which we can find an analytical expression for, by evaluating the leading order term of the slow time tendency equation for the slow mode \(z_0\):

\[\partial_T{z_0} =\boldsymbol p_0 \cdot \boldsymbol S(\boldsymbol z_0, \boldsymbol z_0) \quad \Rightarrow \quad z_{0,0,1} = \boldsymbol p_0 \cdot \boldsymbol I_{0}\]

taking the \(k\)-th derivative yields:

\[z_{0,0,k} = \boldsymbol p_0 \cdot \boldsymbol I_0^{(k-1)}\]

References

[1]

Machenhauer, B. (1977). On the dynamics of gravity oscillations in a shallow water model, with applications to normal mode initialization. Beitr. Phys. Atmos, 50(1).

[2]

Baer, F., & Tribbia, J. J. (1977). On complete filtering of gravity modes through nonlinear initialization. Monthly Weather Review, 105(12), 1536-1539.

[3] (1,2)

Warn, T., Bokhove, O., Shepherd, T.G. and Vallis, G.K. (1995), Rossby number expansions, slaving principles, and balance dynamics. Q.J.R. Meteorol. Soc., 121: 723-739. https://doi.org/10.1002/qj.49712152313

[4]

Eden, C., Chouksey, M., & Olbers, D. (2019). Gravity wave emission by shear instability. Journal of Physical Oceanography, 49(9), 2393-2406.

__init__(mset: ModelSettingsBase, order=3, epsilon=None, use_model=True, time_step_factor=0.01, use_discrete=True, enable_dealiasing=True) → None

Methods

__init__(mset[, order, epsilon, use_model, ...])
bilinear_form(z1, z2)	Calculate the symmetrical bilinear form.
create_state(order[, spectral])	Create a state from the balanced state up to a given order.
interaction(order_series, order_derivative)
reset_fields()	Reset the fields.

create_state(order, spectral=False)

Create a state from the balanced state up to a given order.

Parameters

orderint

The order up to which the state is created.

spectralbool (default: False)

Whether the returned state should be in spectral space.

bilinear_form(z1: VectorField, z2: VectorField) → VectorField

Calculate the symmetrical bilinear form.

Description

The symmetrical bilinear form is defined as:

\[\boldsymbol S = \frac{1}{2} \left( \boldsymbol N(\boldsymbol z_1 + \boldsymbol z_2) - \boldsymbol N(\boldsymbol z_1) - \boldsymbol N(\boldsymbol z_2) \right)\]

where \(\boldsymbol N\) is the nonlinear advection term. The interaction term is a symmetric bilinear form.

Parameters

z1State

The first state (spectral space).

z2State

The second state (spectral space).

Returns

SState

The interaction term (spectral space).

interaction(order_series: int, order_derivative: int) → VectorField

reset_fields()

Reset the fields.