[MITgcm-support] Time dependency of energy for internal gravity waves

[Kevin Ha]

Hi MITgcmers,

I am trying to implement a simple configuration to model some nonlinear
properties of internal gravity waves in the ocean later on. To do so, I am
generating internal waves in a confined box (inviscid ocean with linear EOS
in stratified temperature, no stratified salt) using specific boundary
conditions. However, it seems that the total energy exhibits some spurious
periodical time dependencies that are not supposed to appear according to
the linear theory of internal gravity waves. My guess would be that I might
get a numerical instability somewhere. Here is a detail of my issue.

My reference configuration consists of a 2D x-z stratified non-rotating
domain with nx = 800; ny = 1, nz = 350; dx = dy = 30; dz = 10 (see below
this mail for more details about the configuration extracted from my data
file), with periodic boundary conditions along the x direction and w=0 at
the top and bottom along the z direction, and with the following
simulations parameters:
Nz = 1e-3; %buoyancy frequency
g = 9.81;
alpha = 2e-4; %thermal expansion coefficient
Tref = Nz^2 *(z+Hz/2) /(alpha*g); %stratified temperature
My simulations are run during a five day period (nTimeSteps =  28800,
deltaT = 15.)

The initial setup for the generated wave calls the following definitions
(indicated in the gendata.m file):
nkx0 = 10; %integer for horizontal BC
nkz0 = 20; %integer for vertical BC
kx0 = 2*pi*nkx0/Lx; %periodic BC
kz0 = pi*nkz0/Hz; %imposed by w=0 at top and bottom
k0 = sqrt(kx0^2 + kz0^2); %norm of the wavevector of the generated internal
gravity wave
omega0 = Nz*kx0/k0; %wave frequency

psi0 = sqrt(4*energyfactor)/k0 %wave amplitude with energyfactor = 10^(-4);
U0 = kz0*psi0; %x velocity amplitude
theta0 = - 1/(alpha*g) * (kx0/omega0) * Nz^2 * psi0; %potential temperature
amplitude

In this box, the generated wave has initial condition (indicated in the
gendata.m file):
U = U0.*cos(kz0*z).*cos(kx0*x)
V = zeros(nz,nx);
T = Tref + theta0.*sin(kz0*z).*cos(kx0*x)

The linear theory predicts a resulting stationary wave with time evolution:
U = U0.*cos(kz0*z).*cos(kx0*x-omega0*t)
V = zeros(nz,nx);
W = U0.*kx0./kz0.*sin(kz0*z).*sin(kx0*x-omega0*t)
T = Tref + theta0.*sin(kz0*z).*cos(kx0*x-omega0*t)
W being obtained using the continuity equation.

The mechanical energy of the wave is then given by:
E = KE + PE;
KE = 0.5*(U.^2+V.^2+W.^2) %kinetic energy;
PE = 0.5*(b/N)^2 %potential energy with b the buoyancy linked with the
potential temperature by b= alpha.*g.*theta =
alpha.*g.*theta0.*sin(kz0*z).*cos(kx0*x-omega0*t) (linear EOS);
So all calculations done, the total energy is equal to:
E = 0.5.*psi0.^2.*[kx0.^2.*sin(kz0*z).^2+kz0.^2.*cos(kx0*x-omega0*t)).^2]
Then averaging over the 2D space x-z:
e = <E> = (1/4).*k0.^2.*psi0.^2
Hence the total energy of the wave averaged over space e = <E> should not
depend on time. The same commentary applies for the space mean kinetic and
potential energy ke = <KE> and pe = <PE> and we have equipartition of
kinetic and potential energy as ke = pe = (1/8).*k0.^2.*psi0.^2 (non
rotating domain).

However, my simulations suggest a time dependency of (e,ke,pe), confirmed
when looking at the Fourier transform of the mean space energy e (see
screen capture *EnergyTimeSeries.png* and *EnergyFFT.png*). The
equipartition is also affected even though it is respected in space+time
average (but we are interested in the space average only). Can someone tell
me what could be wrong apparently in my setup? Any ideas on how to correct
that?

You will also find attached the time series of U, V, T and W (
*timeSeries.png*) with their Fourier transform (*FFTtime*), and a video
showing the time evolution of these quantities and of their space Fourier
transform in a (kx,kz) domain (*FFTkxkzrun.avi*)

Thank you in advance,

Kevin

 EnergyFFT.png
<https://drive.google.com/file/d/1UNsA7ut1bbj2-ZA-5rnXzpx7fByfZ3s0/view?usp=drive_web>
 EnergyTimeSeries.png
<https://drive.google.com/file/d/1r5ZvJwehTolIpQkdvWVw9-jJPormdMmJ/view?usp=drive_web>
 FFTkxkzrun.avi
<https://drive.google.com/file/d/1vi9cVrFJWvrp3PibZlpP61AiKxKTqbe3/view?usp=drive_web>
 FFTtime.png
<https://drive.google.com/file/d/13dDtlq0GigcYrqOA4I8HwNaK-l-gtThQ/view?usp=drive_web>
 timeSeries.png
<https://drive.google.com/file/d/1-AuODYPDz24O13AkI6u9rigz8FRb4qoI/view?usp=drive_web>

#########

Config options taken from the data file:

 no_slip_sides  = .TRUE.,
 no_slip_bottom = .TRUE.,

 viscAh  = 0.E-3,
 viscAz  = 0.E-3,

 diffKhT = 0.E-4,
 diffKzT = 0.E-4,

 f0   = 0.E-4,
 beta = 0.E-11,

 tAlpha = 2.E-4,
 sBeta  = 0.E-4,

 gravity = 9.81,
 gBaro   = 9.81,

 rigidLid = .FALSE.,
 implicitFreeSurface=.TRUE.,
 exactConserv = .TRUE.

 nonHydrostatic = .TRUE.,
 useSingleCpuIO = .TRUE.,

 readBinaryPrec  = 64,
 writeBinaryPrec = 64,
 writeStatePrec  = 64,

 staggerTimeStep=.TRUE.

 saltStepping  = .FALSE.,

 tempAdvScheme = 7, !7th Order One Step method with Monotonicity Preserving
Limiter ???? GAD.h

 hFacMin=0.1,

 &

# Elliptic solver parameters
 &PARM02
 cg2dMaxIters       =  10000,
 cg2dTargetResidual = 1.E-14,

 cg3dMaxIters       =    400,
 cg3dTargetResidual = 1.E-14,
 &

# Time stepping parameters
 &PARM03

 niter0 = 0,

 nTimeSteps  =  28800,
 deltaT      =    15.,
 dumpFreq    =  900.,
 monitorFreq = 43200.,

 abEps       =    0.1,
 pChkptFreq  = 86400.,
 chkptFreq   =    0.0,
 &

# Gridding parameters
 &PARM04
 usingCartesianGrid=.TRUE.,
 usingSphericalPolarGrid=.FALSE.,
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.mitgcm.org/pipermail/mitgcm-support/attachments/20210120/d24eb619/attachment-0001.html>