Artificially strong surface velocities (BODYFORCE?)

[aryansafaie]

Dear all,

My grid is a wedge-shaped f-plane domain with a constant bottom slope, initialized with uniform temperature and density fields, and zero u,v,w velocities everywhere.
The only forcing is a constant alongshore surface momentum flux (via ana_smflux), using BODYFORCE to distribute it to all vertical levels (levsfrc = 1), monotonically increasing as a function of cross-shore distance.
The issue is that I am not understanding why the alongshore (u) velocities in the surface cells are so much larger than in the cells just beneath them.

As an example, I've attached a plot showing instantaneous snapshots of depth profiles of alongshore velocity at various cross-shore locations.
The red hues are for a run under downwelling-favorable circulation, and the blue hues are for a run under upwelling-favorable circulation.
Notice the kink at the surface of the velocity profiles; that's the feature that I'm trying to figure out.


I'm not sure where to even begin looking to address this, so I would appreciate any guidance or help.
Thank you very much;


More details are listed below:
dt = 2 s
ndtfast = 35
Coordinates: x positive alongshore (east), y positive cross-shore (north)
theta_s = 3; theta_b = 1; Tcline = 25;
dx = dy = 10 m, bottom slope = 0.03, and the domain is 1 km in the alongshore by 5 km in the cross-shore;
Minimum (onshore) depth = 1 m, max (offshore) depth = 151.3 m

Boundary Conditions:
All onshore (south) boundaries are CLOSED
All alongshore (E/W) boundaries are PERIODIC.
The offshore (north) boundaries are as such:
Che : free-surface
Shc : 2D U-momentum
Shc : 2D V-momentum
Rad : 3D U-momentum
Rad : 3D V-momentum
Gra : mixing TKE
Gra : temperature
Gra : salinity


Maximum grid stiffness ratios: rx0 = 1.034483E-01 (Beckmann and Haidvogel)
rx1 = 5.793344E+00 (Haney)

CPP's:
ANA_BSFLUX Analytical kinematic bottom salinity flux.
ANA_BTFLUX Analytical kinematic bottom temperature flux.
ANA_DRAG_GRID Analytical spatially varying bottom roughness length.
ANA_FSOBC Analytical free-surface boundary conditions.
ANA_INITIAL Analytical initial conditions.
ANA_M2OBC Analytical 2D momentum boundary conditions.
ANA_M3OBC Analytical 3D momentum boundary conditions.
ANA_SMFLUX Analytical kinematic surface momentum flux.
ANA_SSFLUX Analytical kinematic surface salinity flux.
ANA_STFLUX Analytical kinematic surface temperature flux.
ASSUMED_SHAPE Using assumed-shape arrays.
AVERAGES Writing out time-averaged nonlinear model fields.
BODYFORCE Momentum stresses as body-forces.
CANUTO_A Canuto A-stability function formulation.
DIAGNOSTICS_TS Computing and writing tracer diagnostic terms.
DIAGNOSTICS_UV Computing and writing momentum diagnostic terms.
DIFF_GRID Horizontal diffusion coefficient scaled by grid size.
DOUBLE_PRECISION Double precision arithmetic.
GLS_MIXING Generic Length-Scale turbulence closure.
LIMIT_BSTRESS Limit bottom stress to maintain bottom velocity direction.
MIX_GEO_TS Mixing of tracers along geopotential surfaces.
MIX_GEO_UV Mixing of momentum along geopotential surfaces.
MPI MPI distributed-memory configuration.
NONLINEAR Nonlinear Model.
!NONLIN_EOS Linear Equation of State for seawater.
N2S2_HORAVG Horizontal smoothing of buoyancy and shear.
OUT_DOUBLE Double precision output fields in NetCDF files.
PJ_GRADPQ4 Quartic 4 polynomial Pressure Jacobian (Shchepetkin, 2002).
POWER_LAW Power-law shape time-averaging barotropic filter.
PROFILE Time profiling activated .
K_C4ADVECTION Fourth-order centered differences advection of TKE fields.
RADIATION_2D Use tangential phase speed in radiation conditions.
!RST_SINGLE Double precision fields in restart NetCDF file.
SALINITY Using salinity.
SOLVE3D Solving 3D Primitive Equations.
TS_MPDATA Recursive flux corrected MPDATA 3D advection of tracers.
TS_MPDATA_LIMIT Further limiter in upwind corrector fluxes for stability.
TS_DIF2 Harmonic mixing of tracers.
UV_ADV Advection of momentum.
UV_COR Coriolis term.
UV_C4ADVECTION Fourth-order centered differences advection of momentum.
UV_DRAG_GRID Spatially varying bottom roughness length.
UV_LOGDRAG Logarithmic bottom stress.
UV_VIS2 Harmonic mixing of momentum.
VAR_RHO_2D Variable density barotropic mode.
VISC_GRID Horizontal viscosity coefficient scaled by grid size.

[wilkin]

What do you have for LEVBFRC?

[aryansafaie]

Hi John,
1 for levbfrc.

[wilkin]

It's not obvious to me why you would be getting that surface shear effect, and quickly looking at the BODYFORCE code I don't see any problems with the k indexing.

Maybe try stripping out a lot of the options you have that you don't need.

In particular, MIX_GEO_UV (and MIX_GEO_TS) - you have no density gradients so no worries about diapycnal mixing.

Your grid resolution is uniform so you don't need DIFF_GRID or VISC_GRID.

PJ_GRADPQ4 is an unconventional choice of pressure gradient - is there a reason you picked this?

The term that might be monkeying with the surface is the GLS_MIXING. Perhaps try using ANA_VMIX instead of a closure scheme. With BODYFORCE you don't need the a closure scheme to mix the wind stress momentum into the water column. It's only for subsequent shear in the flow.

[aryansafaie]

Hi John,

Thanks for the advice.
Now this is interesting: I was inspired by what you mentioned regarding my use of the k-epsilon closure; so I ran another quick simulation, this time using KPP instead.
A comparison of vertical profiles of alongshore velocities and eddy viscosity at various cross-shore locations is shown below.
Unlike the GLS run, the KPP run lacks the high surface shear, and at the surface, the eddy viscosity goes exactly to 0 (unlike k-epsilon).
Could part of the issue be that I have the Charnok CPP turned off, even though I'm prescribing a surface "wind stress" via ana_smflux?
I'm going to go down that path and I'll report back.

Thanks again!
-aryan

[aryansafaie]

I believe it's about time I post an update:

After reexamining Warner et al 2005 (OM), it seemed like a good idea to reconsider the surface boundary condition for TKE (equation 51 in that paper).
The file gls_corstep.F contains code to implement the surface BC for TKE; from my prior simulations, it seems that the surface TKE boundary condition was too large: as per k_s = u_*²/(c_mu⁰)², it seems that all of the surface momentum flux I'm prescribing via ana_smflux.h is going into that surface cell, and not just the scaled u_*surface² that I would expect via the BODYFORCE CPP.

Therefore, I ran two additional simulations in which I manipulate gls_corstep:
#1: I set TKE@the surface = 0
#2: I multiply the surface TKE boundary condition by the ratio of the grid cell thickness to the water column depth:
TKE_s = (∆z/H_total)*[(u_*)/(c_mu⁰)]²
this latter idea is to ensure that the TKE surface boundary condition takes into account the BODYFORCE CPP.

I've attached plots showing vertical profiles at various cross-shore locations of eddy viscosity, alongshore velocity, and TKE; the black line is the default GLS k-epsilon scheme, the red line is just setting TKE_surface = 0, the green line is my little hack (#2 described above), and the blue line is using the KPP scheme as a sort of sanity check.
It seems that in either of my little manipulations to gls_corstep.F, more physically realistic results are obtained by (correctly?) reducing the TKE at the surface.

These results are fine and dandy for my purposes, but I just wanted to be sure I'm not overlooking anything or making a horrible mistake.

Thanks again,
-aryan