Mean Flow from Phase Averages in the 2D Boussinesq Equations

Abstract

The atmosphere and ocean are described by highly oscillatory PDEs that challenge both our understanding of their dynamics and their numerical approximation. This paper presents a preliminary numerical study of one type of phase averaging applied to mean flows in the 2D Boussinesq equations that also has application to numerical methods. The phase averaging technique, well-known in dynamical systems theory, relies on a mapping using the exponential operator, and then an averaging over the phase. The exponential operator has connections to the Craya–Herring basis pioneered by Jack Herring to study the fluid dynamics of oscillatory, nonlinear fluid dynamics. In this paper, we perform numerical experiments to study the effect of this averaging technique on the time evolution of the solution. We explore its potential as a definition for mean flows. We also show that, as expected from theory, the phase-averaging method can reduce the magnitude of the time rate of change in the PDEs, making them potentially suitable for time stepping methods.

1. Introduction

The mean flow in the present study is a time-average of a system of equations that has been mapped using the exponential operator. The idea is that oscillatory dynamics, once they are mapped into a new ‘moving frame’, can be composed of a slowly varying mean and fast departures from that mean. This idea has a long history in studying the time evolution of the dynamics and in the numerical approximation of oscillatory ODEs and PDEs, which we discuss in the next paragraph. The technique explored in this paper has also been proposed as a way of time stepping that uses a serial prediction step and a time-parallel correction step. The time stepping technique shows promise but has only been applied to ODEs and 1D PDEs. In this paper, we use the mapping and averaging technique as a diagnostic in numerical experiments to (1) study the impact of this mapping and averaging on the time-evolution of the 2D Boussinesq equations, (2) to provide an initial assessment of whether it might be a useful definition of a mean flow, and (3) to assess its potential usefulness in predictor–corrector-type numerical methods.

The technique of mapping and averaging, now a classic method in dynamical systems analysis, is discussed and developed in [1], where they also give a brief history of the averaging method since the lectures of Jacobi, which date to 1842. One of their key points is that the oscillatory system of equations needs to be mapped into their ‘standard form’. This mapping is equivalent to our use of the exponential map in this paper. Further, the eigenbasis used for the exponential map in this paper has a relationship to the Craya–Herring basis pioneered by Jack Herring in [2]. Furthermore, in several papers, including [3,4], Herring and his collaborators used this basis as a means to understand the dynamics of waves relative to vortical aspects of the flow. In this paper we explore decompositions of this nature, but for constructing phase averages of the nonlinear dynamics. Ref. [1] gives a demonstration (page 28) of why straightforward averaging without mapping, which they call naive averaging, gives approximation errors. The technique in this paper also plays a role in the theory of fast singular limits [5,6,7,8,9,10]. The mapping concept without averaging has also been useful in the study of geophysical fluid dynamics, e.g., [11,12,13], to study nonlinear resonances.

Another way of understanding this technique, for equations in the form of (1), is to think of it as nonlinear phase averaging because, once the PDE has been mapped using the exponential operator (2), all the phase information arises in the nonlinear term (see Equation (3)). Recent papers that explicitly discuss phase averaging in geophysical fluid dynamics include [14,15,16].

From a PDEs point of view, the averaging technique in this paper is related to the concept of ‘cancellation of oscillations’ [7,10] to treat hyperbolic PDEs. The authors of [8] refer to this technique and treat geophysical applications, that is, the Rotating Boussinesq Equations and the Rotating Shallow Water Equations. In this context, an exponential operator and infinite time averaging leads to limiting equations for the ‘slow’dynamics. This type of averaging was further developed in [17,18] for geophysical applications. In [18,19], the analysis was extended to the case when more than one fast time scale is present.

The present study, like those mentioned in the previous paragraph, uses an exponential operator as a mapping. We emphasize that using the exponential operator mitigates the linear oscillations in the problem, but the oscillations are not eliminated—there are oscillations that develop from the nonlinear terms. Therefore, as an additional ingredient, we evaluate an integral with the scaled kernel function $\mathcal{K}$, that is, we compute a convolution. It is this technique that potentially can describe general time-mean flows and that can potentially be used in numerical time-stepping methods. It differs from the technique used in the study of fast singular limits in that the averaging is over a finite window and uses a smooth integration kernel over just a few oscillations, rather than an infinite time average.

Numerically, the classical method of averaging is described in [1]. In geophysical fluid dynamics applications, a pioneering paper [20] directly applies the method of [5,6] to the shallow water equations; they found that the method worked well when the oscillations were fast, but became less accurate as the frequency of the oscillations decreased. The phase-averaged mean flows studied in this paper use a technique proposed by [21,22,23] as a coarse time stepping algorithm for a numerical time-stepping method. This technique has also been used for time-stepping on the sphere [16].

The oscillatory equations we study in this paper are the 2D Boussinesq equations given in (6)–(8); before we discuss those, we first describe the phase-averaging method we use in this paper. We therefore study the more general form of oscillatory equations:

$$\frac{d\mathbf{u}}{dt}+\frac{1}{\epsilon }\mathcal{L}\mathbf{u}=\mathcal{N}\left(\mathbf{u}\right).$$

(1)

The linear operator $\mathcal{L}$ has purely imaginary eigenvalues and makes the problems stiff. Additionally, the non-linearity $\mathcal{N}$ is bi-linear in fluid applications. One possibility is to eliminate the linear term in the time evolution problem (1) by using the transformation

$$\mathbf{w}\left(t\right)=exp\left(\frac{\mathcal{L}}{\epsilon }t\right)\mathbf{u}\left(t\right).$$

(2)

This leads to the following problem:

$$\frac{d\mathbf{w}}{dt}=exp\left(\frac{\mathcal{L}}{\epsilon }t\right)\mathcal{N}\left(exp\left(-\frac{\mathcal{L}}{\epsilon }t\right)\mathbf{w}\right).$$

(3)

Applying averaging to further mitigate the oscillatory stiffness leads to a new problem:

$$\frac{d\tilde{\mathbf{w}}}{dt}=\frac{1}{\eta }{\int }_{-\eta /2}^{\eta /2}\mathcal{K}\left(\frac{s}{\eta }\right)exp\left(\frac{\mathcal{L}}{ϵ}(s+t)\right)\mathcal{N}\left(exp\left(-\frac{\mathcal{L}}{ϵ}(s+t)\right)\tilde{\mathbf{w}}\left(t\right)\right)ds.$$

(4)

Solving Equation (4) can be very useful for numerical time-stepping, such as in the Parareal method, because the time rate-of-change of the unknowns becomes less oscillatory. This path was extensively studied in [21,22,23].

From a fluid dynamics point of view, this type of equations may be interpreted as the mean flow. In the present study, we want to step in a new, less-explored direction by applying the transformation and averaging to a given solution $\mathbf{u}\left(t\right)$ of problem (1) as a diagnostic tool and investigating if such formulations provide a suitable formulation for mean flows. In this case, we evaluate

$$\overline{\mathbf{w}}\left(t\right)=\frac{1}{\eta }{\int }_{-\eta /2}^{\eta /2}\mathcal{K}\left(\frac{s}{\eta }\right)\mathbf{w}(t+s)ds,$$

(5)

where $\mathbf{w}$ is given by Relation (2). The difference between Equations (4) and (5) is that, in the first equation, the averaging procedure is applied to the right-hand side of a differential equation, whereas, in the second equation, it is applied to a given solution of that differential equation. Thus, in the first case, averaging is applied to make a time evolution problem more well-behaved for numerical computations; in the second case, it serves to extract information from data. The averaging is discussed in more detail in Section 2.3.

Already, in [8], it was suggested to use exponential mapping for the formulation of the mean flow for geophysical fluid dynamics applications. However, the approach described there differs from the ideas in the present expositions in the following ways. First, the exponential in this exposition is applied to data representing a solution to the Boussinesq equations. We do not use it to derive analytical expressions for leading order solutions when the method of multiple scales is applied. Second, for the formulation of the mean flow, the authors of [8] proposed to consider only modes which satisfy the condition for resonant triads. In the present exposition, this constraint is relaxed and not only resonant modes, but also other slow modes are admitted for the formulation of the mean flow.

In the present study, the focus is the 2D Boussinesq equations. In particular, we study the case of an initial buoyancy distribution that decays with time:

The Boussinesq equations are

$$\begin{array}{ccc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}{u}_t+u{u}_x+w{u}_z+\frac{\partial }{\partial x}\left({\Delta }^{-1}\left(-\nabla ·\left(\mathit{u}·\nabla \mathit{u}\right)-N\frac{\partial {\rho }^{\prime }}{\partial z}\right)\right)& =& \nu \Delta u,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill w_t+u{w}_x+w{w}_z+\frac{\partial }{\partial z}\left({\Delta }^{-1}\left(-\nabla ·\left(\mathit{u}·\nabla \mathit{u}\right)-N\frac{\partial {\rho }^{\prime }}{\partial z}\right)\right)+N{\rho }^{\prime }& =& \nu \Delta w,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{\partial {\rho }^{\prime }}{\partial t}+u{\rho }_x^{\prime }+w{\rho }_z^{\prime }-Nw& =& \kappa \Delta {\rho }^{\prime },\hfill \end{array}$$

(8)

where the the Brunt–Väisälä frequency is $N=\sqrt{\frac{bg}{{\rho }_0}}$. We assume periodic boundary conditions with the vector of unknowns $\mathbf{u}=(\mathbf{v},{\rho }^{\prime })$. The divergence constraint is given by $div\mathbf{v}=0$, where $\mathbf{v}={(u,w)}^T$. The pressure obeys a Poisson equation (see [24]). We study the averaging technique for oscillating and decaying dynamics.

In [8], the following identity was derived for a leading order solution of the rotating Boussinesq equations:

$${\mathbf{u}}^0=e^{-\tau \mathcal{L}}\overline{\mathbf{u}},$$

(9)

where $\overline{\mathbf{u}}$ satisfies an averaged equation, which contains the non-linearity and is given by

$$\frac{d\mathbf{u}}{dt}=-\underset{\tau \to \infty }{lim}\frac{1}{\tau }{\int }_0^{\tau }{e}^{s\mathcal{L}}\mathcal{B}(e^{s\mathcal{L}}\mathbf{u}\left(t\right),e^{s\mathcal{L}}\mathbf{u}\left(t\right))ds.$$

(10)

Note that we take $\tau \to \infty$ in Equation (10). The integrand in Equation (10) can be represented using a Fourier decomposition in space. Applying the Riemann–Lebesque lemma, one sees that only the modes which form direct resonances are retained by the intregration process. However, we relax this constraint and allow some slow oscillations too. This provides a new tool to define mean flows, where we allow selected oscillations to be in the formulation. The new averaging procedure, which differs from the one given by Equation (10), is of central significance to decide which oscillations remain in the solution and is explained in more detail in Section 2.3.

In the remainder of the article, we describe the details of the diagnostic we use to compute the time-mean flow. In the results section, we (1) present preliminary numerical results that show the behavior of the mapping and averaging on the time-evolution of the buoyancy, (2) characterize the mean flow that results from the mapping and averaging, including exploring the effect of changing the averaging window, and (3) compute the time rate-of-change of the buoyancy, to assess whether it can mitigate the oscillatory stiffness for use in numerical methods.

2. Methods

We compute a temporal mean flow for solutions to the Boussinesq Equations. For that, the following procedure is pursued. First, a solution to the Boussinesq Equations, denoted as $\mathbf{u}$, is computed. Particulars can be found in Section 2.1. Then, the exponential of the linear operator is applied to the solution and we obtain the mapped solution $\mathbf{w}$. The computation of the exponential is described in Section 2.2 and Appendix A. In the next step, the mapped solution is averaged according to Equation (5). Thus, we compute $\overline{\mathbf{w}}$. The mapped solution $\mathbf{w}$ can be represented in terms of spatial Fourier modes. When applying the averaging from Equation (5), we can choose, which modes shall be averaged or remain in the solution, by adjusting $\eta$, denoted as the averaging window. If the period of a mode is larger than $\eta$, the mode will still be present in the solution—see Section 2.3 for details. In the last step, we apply the inverse of the transformation (2).

Performing the previously described procedure, we create several data sets:

• The unchanged data from the numerical solution of the 2D Boussinesq equations $\mathbf{u}$;
• The data after applying the transformation, yielding the mapped solution $\mathbf{w}$;
• The data after applying the transformation and the averaging procedure, which give the averaged solution $\overline{\mathbf{w}}$;
• The data after applying the inverse transformation $\overline{\mathbf{u}}$.

The above items are illustrated in Figure 1. Comparison of the data sets helps to assess how well the proposed procedure is suited to describe mean flows and potentially to be used for time-stepping methods.

2.1. Numerical Computation of Solutions to the Boussinesq Equations

We compute numerical solutions using the Dedalus software [25]. We use a Fourier discretization in space in a $2\pi$ periodic 2D domain and 3rd-order 4-stage Runge–Kutta Method (RK433) time-stepper. For all the simulations, the time step is $\Delta$t = 0.001, the spatial grid is $(N_x,N_z)$ with $Nx=Nz=128$, the viscosity is $\nu$ = 0.0001, and the dissipation for buoyancy is $\kappa$ = 0.0001. For the initial condition, $u(0,\mathbf{x})=w(0,\mathbf{x})=0$, while the initial buoyancy is

$${\rho }^{\prime }(0,\mathbf{x})=\sqrt{{ϵ}_f}\sin \left(3z\right)\sin \left(3x\right),$$

(11)

with ${ϵ}_f=25$. The energy oscillates between kinetic and potential energy in time while it decays. The higher the value of N, the more frequent the oscillations. Figure 2 shows an example of the dynamics in decay as it oscillates between kinetic and potential energy. At some time after the simulation begins, higher frequency oscillations are excited, providing a challenge for studying oscillatory mean flows.

In this study, we use 3 different values of N: $N=1,10,$ and 20. We define the Froude number to be a function of the character of the initial condition, using similar ideas to [12], who studied the forced case

$$Fr=2\frac{{\left({ϵ}_f{k}_f^2\right)}^{\frac{1}{3}}}{N},$$

(12)

which gives us equivalent Froude numbers of $Fr=10.6,1.06,$ and 0.53.

2.2. Fourier Analysis and the Computation of the Exponential of $\mathcal{L}$ with the Divergence Constraint

For the implementation of the exponential of the linear operator, we apply a Fourier expansion in space to the solutions of the Boussinesq equations, i.e., we have

$$\mathbf{u}(t,\mathbf{x})=\sum _{k=-N}^N\widehat{\mathbf{u}}\left(t\right)e^{i\mathbf{k}·\mathbf{x}}.$$

(13)

Each mode $e^{i\mathbf{k}·\mathbf{x}}$ is applied to the linear operator, which is given by

$$L\mathbf{u}=\left(\begin{array}{c}-\frac{\partial }{\partial x}\left({\Delta }^{-1}N\frac{\partial {\rho }^{\prime }}{\partial z}\right)\\ -\frac{\partial }{\partial z}\left({\Delta }^{-1}N{\rho }^{\prime }\right)+N{\rho }^{\prime }\\ -Nw.\end{array}\right)$$

(14)

Consequently, a $3\times 3$ matrix is obtained for each mode. Then, the matrix exponentials for the different modes are computed. To compute the exponential of the linear operator, it is sufficient to know the exponential for each mode $e^{i\mathbf{k}·\mathbf{x}}$. Due to linearity, we have

$$e^{\frac{\mathcal{L}}{ϵ}t}u(t,\mathbf{x})=\sum _{k=-N}^N{e}^{\frac{\mathcal{L}}{ϵ}t}{e}^{i\mathbf{k}·\mathbf{x}}\widehat{\mathbf{u}}\left(t\right)=\sum _{k=-N}^N{B}_k\widehat{\mathbf{u}}\left(t\right)e^{i\mathbf{k}·\mathbf{x}},$$

(15)

where $B_k{e}^{i\mathbf{k}·\mathbf{x}}=e^{\frac{\mathcal{L}}{ϵ}t}{e}^{i\mathbf{k}·\mathbf{x}}$. Applying the exponential, we obtain $\mathbf{w}$, given by (2).

Here, we briefly describe the computation of the matrix exponential. Further details of the computations can be found in Appendix A. The horizontal wave number is denoted as $k_1$. With $k_3$, we refer to the vertical wave number. First, the eigenvalues are computed. Then, the corresponding eigenvectors are identified. With the eigenvectors, we know the transformation matrix $T_k$ and can also compute its inverse $T_k^{-1}$. In the computations, we distinguish between four different cases (first: $k_1=0,k_3=0$; second: $k_1=0,k_3\ne 0$; third: $k_1\ne 0,k_3=0$; fourth: $k_1\ne 0,k_3\ne 0$). The matrices are diagonalizable, except for the case $k_1=0,k_3\ne 0$. For that case, we compute the Jordan normal form. To compute the exponential of the matrix $B_k$, where $B_k$ comes from the spatial discretization of the linear operator and depends on the wave numbers, we use the identity $exp\left(B_k\right)=T_{k}exp\left(M_k\right)T_k^{-1}$. In the relation, the matrix $M_k$ is the diagonal matrix or the matrix in Jordan normal form with the computed eigenvalues. Finally, we enforce the divergence constraint.

2.3. Averaging with a Convolution Integral

A transformed solution of the Boussinesq equations $\mathbf{w}$ is convolved by applying a scaled filter function $\mathcal{K}$:

$$\mathcal{K}\left(\frac{s}{\eta }\right)=\frac{1}{K_0}exp\left(\frac{1}{(s/\eta -1/2)(s/\eta +1/2)}\right),$$

(16)

where the constant $K_0$ is chosen so that integrating the kernal function along the real line gives 1. We scale the filter function by choosing an averaging window $\eta$. The idea behind the choice of the averaging window $\eta$ is that it must mitigate the fast oscillations while leaving the coarser behavior of the solution unaffected. The bigger the averaging window, the more oscillations are mitigated. The smaller the averaging window, the more oscillations are resolved. Temporal averaging was also investigated in other contexts, like ordinary differential equations [1], heterogeneous multiscale methods [26,27], and PDE analysis [8].

To further illustrate the averaging, let us suppose that a slow function is superimposed by a fast periodic function with zero mean. When we compute the integral of the fast function over an interval of length $\eta$, assuming that $\eta$ is as large as a few times the period of the fast periodic function, the positive and negative contributions cancel each other. However, the integrand is weighted by the scaled kernel function $\mathcal{K}$ with compact support, which decays fast close to the boundary of the compact support. In addition, we do not assume that the length of the exact period is known. Thus, in the general case, the oscillations are not exactly canceled, but rather they are mitigated, as can be seen in Figure 3. Technical details of the averaging can be found in [27] in Lemma 2.2.

3. Results and Discussion

We begin with a description of the time evolution of the dynamics (unaveraged) for the specific case studied in this work. In all cases, the initial condition is given by (11) and the velocity components are zero initially.

Figure 2 shows the evolution of total, potential, and kinetic energy for $N=10$ on the left and $N=1$ on the right. The data are shown in black lines for the unaveraged case with total energy: solid, kinetic energy: dotted, and potential energy: dashed. These graphs show that the dynamics decay as it evolves, and that the case with the larger Brunt–Väisälä frequency has more frequent energy exchanges between potential and kinetic energy compared to the $N=1$ case.

These more frequent energy exchanges in potential and kinetic energy can also be seen by examining the time history of $\mathbf{u}\left[3\right](x_p,z_p,t)={\rho }^{\prime }(x_p,z_p,t)$ at a single, randomly chosen point $(x_p,z_p)=(\frac{2\pi 60}{128},\frac{2\pi 60}{128})$. Other choices of the location of the point show similar results. This will be explored in the next section where we examine vertical slices, which show the time evolution for all the values of z at a fixed $x_p$. Compare the solid black line in the upper right panel of Figure 4 for $N=10$ to the same style of line in Figure 4 for the $N=1$ data. For the $N=10$ case, the solution has more frequent oscillations than the $N=1$ case. For both values of N, the solution breaks up into higher oscillations as the solution decays. This feature is used to illustrate the effect of the averaging in the following sections.

Finally, before we go on to study the effect of the averaging, we look at the time evolution of the single point data in the mapped domain $\mathbf{w}\left[3\right](x_p,z_p,t)$. These data are computed by applying the matrix exponential, as in (2), to the vector $\mathbf{u}$. This has the effect of showing what the time evolution of the solution looks like in the moving frame and allows us to observe the effect of the ‘phase scrambling’ in the nonlinear term, as in (3). The evolution in the moving frame is shown in the solid black line in the left panels of Figure 4. For $N=20$ and $N=10$, there are more frequent oscillations in the moving frame than in the ordinary frame shown in the right panel. In contrast, for the $N=1$ case, the difference between the solution in the moving frame $\mathbf{w}$ and the ordinary frame $\mathbf{u}$ is not as different as the other cases, which is because the moving frame is slow for $N=1$.

In the next three subsections, we examine the effect of the mapping and averaging on these oscillatory, decaying solutions.

3.1. Impact of Mapping and Averaging

We now describe the impact of the mapping and averaging technique on the time evolution of the ${\rho }^{\prime }$ dynamics by examining time series for three different values of N in the context of the four data sets ($\mathbf{u},\mathbf{w},\overline{\mathbf{w}},\overline{\mathbf{u}}$) described in the second paragraph of Section 2. Because we focus on ${\rho }^{\prime }$, when $\mathbf{u},\mathbf{w},\overline{\mathbf{w}},\overline{\mathbf{u}}$ appear in a figure, these are the third component of the vector field.

Figure 4 presents the time series at a single point, with each row of the figure corresponding to a different value of N. In all three rows, the left panel shows the time evolution in the moving frame $\mathbf{w}$ and its average $\overline{\mathbf{w}}$. In contrast, the right panel shows the ordinary variable $\mathbf{u}$ and the average in the moving frame mapped back to the ordinary frame $\overline{\mathbf{u}}$. The significance of the left panels is that they show the result of how the oscillations are being ‘scrambled’ and evolved by the nonlinear term in (2) as the solution undergoes decay from its initial condition. A first observation is that, for the two higher values of N studied, the oscillations in $\mathbf{w}$, the moving frame variables, are faster than those observed in $\mathbf{u}$. For some parameter regimes in 2D Boussinesq dynamics, one might expect the dynamics of $\mathbf{w}$ to be of lower frequency than $\mathbf{u}$, not higher, so further studies of the case with forcing, where low frequency buoyancy layers are formed, could be interesting (see [12]).

A second observation is that, for the higher values of N at later times, when the faster oscillations appear, the averaging from the nonlinear phase scrambler in the moving frame, $\overline{\mathbf{w}}$, tracks the mean value, as expected. This is in contrast to the panel on the right where, for example, for $N=20$ in the range $3.5<t<4$, when the average is mapped back to the ordinary frame, the average does not track the mean value. This is also as expected because, in this domain, the average is modulated by the oscillations. As mentioned in the introduction, ref. [1] gives an example of why straightforward time averages, in the unmapped equations, do not give accurate results. The comparison that we just described, where the modulated mean, $\overline{\mathbf{u}}$, does not necessarily track the mean of $\mathbf{u}$, is a characteristic feature of averages taken in the modulated domain. We discuss this more in the next section where we consider the mean-flow possibilities.

Finally, we make a few observations about the effect of the averaging window $\eta$ on the mean of the time series. In Figure 4, on the left panels, the effect of varying the averaging window is shown. The longer the averaging window, the more the amplitudes are damped. At early times, for $N=20$, the main effect is to reduce the magnitude of the oscillations, but, at later times, the red line, for which a small averaging window $\eta$ = 0.05 was used, shows little difference from the unaveraged case. However, when the largest averaging window is applied, for $\eta$ = 0.2, the blue line for $3<t<4$ shows that the mapping and averaging tracks the mean of the oscillations. This similar behavior can be seen in the unmapped domain on the right, where, as pointed out earlier, the blue line, with the larger averaging window, is slower, but does not follow the mean of the oscillations. For the middle panel, when $N=10$, we show three different averaging window lengths ($\eta$ = 0.05, 0.1, and 0.2). For the largest averaging window (blue line), the mean tracks the mean of the oscillations. This is true in the right panel as well. Finally, for the lowest value of N, $N=1$, the choice of averaging window has a similar effect, except that, at early times, it does not have that much of an effect, while, at later times, when the faster oscillations begin, it appears to track through the mean value.

Figure 5 and Figure 6 show the time evolution of four data sets for every point in the z plane for a fixed $x=x_p$, for $N=10$ and $N=1$, respectively. Rather than grouping $\mathbf{w}$ and $\overline{\mathbf{w}}$ together and $\mathbf{u}$ and $\overline{\mathbf{u}}$ together, as in the previous paragraphs, this time we show a panel for each data set in the same counterclockwise order as shown in Figure 1. The upper left graph shows the time evolution of $\mathbf{u}$, the lower left shows the time evolution in the moving frame $\mathbf{w}$, and the lower right panel shows the average in the moving frame $\overline{\mathbf{w}}$. For $N=10$, $\mathbf{u}$ oscillates between positive and negative values of the buoyancy and then we see the creation of higher frequency waves as the solution decays. We observe that, even when the higher frequency oscillations appear, larger scale (in time) patterns can also be seen. The lower right panel shows the effect of averaging in the moving frame, where the larger structures in time now look smoother. For the case when $N=1$, the oscillatory pattern between positive and negative values are less frequent, since N is smaller. Comparing the two cases $N=10$ and $N=1$, the graphs show that $N=1$ has an overall large-scale (in time) pattern of oscillations, whereas the $N=10$ case has frequent temporal oscillations.

Finally, we examine the effect of the mapping and averaging on the energy. Figure 2 shows the effect of the averaging on the kinetic, potential, and total energy for the $\mathbf{u}$ and $\overline{\mathbf{u}}$ data sets for $N=10$ on the left and $N=1$ on the right. For $N=10$, the mean tracks the oscillation between kinetic and potential, and shows an oscillation in the overall mean flow, which is discussed a bit more in the next section. For $N=1$, the mean total energy decays faster than total energy, due to the average over the oscillations.

3.2. Definition of a Mean Flow

For a preliminary assessment of whether the mapping and averaging procedure could be a suitable definition of a mean flow for the 2D Boussinesq equations, we break the unknowns into its mean and fluctuation:

$$\mathbf{u}=\overline{\mathbf{u}}+{\mathbf{u}}^{\prime },$$

(17)

and in the moving frame:

$$\mathbf{w}=\overline{\mathbf{w}}+{\mathbf{w}}^{\prime }.$$

(18)

Integration over the oscillations would give us zero as positive and negative components cancel each other. From a mean flow, we would expect that only the mean behavior without the departures is retained.

Revisiting Figure 4, we make the following two observations: (1) The fastest variations in the solution are not visible anymore; (2) The amplitudes are damped when the solutions have steep gradients. Therefore, the component which is not retained contains fast oscillations and has positive and negative components.

Figure 7 shows the point-wise departure from the mean in the moving frame, ${\mathbf{w}}^{\prime }$, for three different averaging windows. The main variations from the mean occurs when the filtering window is larger, as we expect, and also where there are higher oscillations at later times in the simulation. We also observe that, once the higher frequency oscillations have begun, the variation from the mean is also nearly periodic. Variation from the mean for the ordinary frame ${\mathbf{u}}^{\prime }$ (not shown) shows a different detailed pattern with a similar periodicity.

Additionally, the choice of the averaging window allows us to determine which oscillations should be filtered and how strong the damping should be. For large averaging windows, more oscillations are filtered and the damping is stronger than for small averaging windows. This means that the definitions of the mean flow and the departures from the mean depend on the choice of the length of the averaging window. In both the moving frame and the ordinary frame, periodic patterns can be seen in the deviation from the mean.

Finally, we point out that total energy in Figure 2 for $N=10$ shows an oscillation in the mean and that the the distance from the mean to the instantaneous energy decay is the energy from the variations from the mean. This oscillatory pattern in the mean flow is a potential area to explore in further work.

Our observation in the last section, that the mapping and averaging tracks the mean in the moving frame but not in the ordinary frame, will also be interesting in the context of fast singular limits, where it is expected that low-frequency structures will appear over time, rather than only high frequencies. To obtain a more complete understanding a wider range of behaviors, such as the creation of low frequencies, would be useful.

3.3. Potential for Use in Numerical Methods

As mentioned in the introduction, the technique explored in this paper has been proposed for use as a prediction step for time-parallel predictor–corrector methods used to solve PDEs of type (1). An error bound for a predictor-step algorithm using this method can be found in [22] and does not require that $ϵ$ be small. As a part of the error estimate, the error bound for the approximation of the time rate of change of the unknowns depends on the averaging window, $\eta$. We therefore are interested in whether the time rate of change of the variable examined in this paper has this property. To study this, we apply a first-order finite difference method to the moving frame $\mathbf{w}$ and $\overline{\mathbf{w}}$ data sets. Figure 8 shows the time rate of change for the $N=10$ case depicted in the time series of Figure 4. This shows that the nonlinear averaging process reduces the amplitude and frequency of the oscillations, making it possible that this method could provide a good predictor for numerical methods. However, it is clear that the filter width matters, in line with the error estimate, and that it may need to be dynamically adjusted. For example, the solid blue line in Figure 8 shows that it has fewer large oscillations than the unaveraged case, but whether this is adequate to take larger time steps is unclear.

We have also pointed out, in previous sections, that, for the simulations in this study, the actual mean tracks $\overline{\mathbf{w}}$ rather than $\overline{\mathbf{u}}$. It would be of interest to explore this more fully in future work.

4. Conclusions

In this work, we have presented numerical studies of the oscillating and decaying 2D Bousinnesq equations and developed a diagnostics to study a nonlinear averaging procedure. The goal was to present a preliminary study of the impact of the mapping and averaging technique by using it as a diagnostic, then assess whether it could make a good definition of a mean flow and potentially help mitigate oscillatory stiffness in numerical methods.

We made the following observations:

• The nonlinear averaging procedure depends strongly on the value of $\eta$, but could potentially provide a reasonable representation of the time-mean of the solutions, whether or not N is large. Further, the average ‘modulates’ the waves, as in (2). This has the effect that the mean $\overline{\mathbf{u}}$ depends on how fast the oscillations are in the exponential operator (2), which could be interesting for further studies. Further work is needed to assess whether this would be a good definition of a mean flow in fluid dynamics. For example, it would be interesting to examine a case that admits a low-frequency solution such as the 2D Boussinesq equations with forcing, as in [12] and examples shown in [8];
• We also observed that the method could potentially be used to take larger time steps, as has been shown in simple examples [21], and that the ability to do so critically depends on the value of $\eta$.