Problems with unforced model -- unrealistic circulation

[Arvin_Yin]

Hi,all
I am modeling the M2 internal tide in South China Sea. First I run the model without any forcing in the domain. The initial T/S are horizontally uniform. Initial velocity and ssh are set to zero. Then I run the model for two weeks. I thought that the ocean would remains at rest when unforced, but finally the model come out a circulation (see the following picture).

When I force the model with M2 tide, the suface Eulerian residual current come out a similar circulation. So, I think it does matters.
I wonder if there are something wrong with the model configuration. And my model configuration is:

Activated C-preprocessing Options:

TIDAL_FORCING
ANA_BSFLUX
ANA_BTFLUX
ANA_FSOBC
ANA_M2OBC
ANA_SMFLUX
ANA_SSFLUX
ANA_STFLUX
ASSUMED_SHAPE
DIFF_GRID
DJ_GRADPS
DOUBLE_PRECISION
GLS_MIXING
KANTHA_CLAYSON
MASKING
MIX_GEO_TS
MIX_S_UV
MPI
NONLINEAR
NONLIN_EOS
N2S2_HORAVG
POWER_LAW
PROFILE
K_GSCHEME
!RST_SINGLE
SALINITY
SOLVE3D
SPLINES
TS_U3HADVECTION
TS_C4VADVECTION
TS_DIF2
UV_ADV
UV_COR
UV_U3HADVECTION
UV_C4VADVECTION
UV_QDRAG
UV_VIS2
VAR_RHO_2D
VISC_GRID

Variable Grid West Edge South Edge East Edge North Edge
--------- ---- ---------- ---------- ---------- ----------

zeta 1 Closed Chapman Exp Chapman Exp Chapman Exp

ubar 1 Closed Shchepetkin Shchepetkin Shchepetkin

vbar 1 Closed Shchepetkin Shchepetkin Shchepetkin

u 1 Closed Radiation Radiation Radiation

v 1 Closed Radiation Radiation Radiation

temp 1 Closed Radiation Radiation Radiation

salt 1 Closed Radiation Radiation Radiation

tke 1 Closed Closed Closed Closed

40320 ntimes Number of timesteps for 3-D equations.
30.000 dt Timestep size (s) for 3-D equations.
20 ndtfast Number of timesteps for 2-D equations between
each 3D timestep.

1.0000E-01 nl_tnu2(01) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 01: temp
1.0000E-01 nl_tnu2(02) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 02: salt
1.0000E+00 nl_visc2 NLM Horizontal, harmonic mixing coefficient
(m2/s) for momentum.

1.0000E-06 Akt_bak(01) Background vertical mixing coefficient (m2/s)
for tracer 01: temp
1.0000E-06 Akt_bak(02) Background vertical mixing coefficient (m2/s)
for tracer 02: salt
1.0000E-05 Akv_bak Background vertical mixing coefficient (m2/s)
for momentum.
5.0000E-06 Akk_bak Background vertical mixing coefficient (m2/s)
for turbulent energy.
5.0000E-06 Akp_bak Background vertical mixing coefficient (m2/s)
for turbulent generic statistical field.

3.000 gls_p GLS stability exponent.
1.500 gls_m GLS turbulent kinetic energy exponent.
-1.000 gls_n GLS turbulent length scale exponent.
7.6000E-06 gls_Kmin GLS minimum value of turbulent kinetic energy.
1.0000E-12 gls_Pmin GLS minimum value of dissipation.
5.4770E-01 gls_cmu0 GLS stability coefficient.
1.4400E+00 gls_c1 GLS shear production coefficient.
1.9200E+00 gls_c2 GLS dissipation coefficient.
-4.0000E-01 gls_c3m GLS stable buoyancy production coefficient.
1.0000E+00 gls_c3p GLS unstable buoyancy production coefficient.
1.0000E+00 gls_sigk GLS constant Schmidt number for TKE.
1.3000E+00 gls_sigp GLS constant Schmidt number for PSI.

1400.000 charnok_alpha Charnok factor for Zos calculation.
0.500 zos_hsig_alpha Factor for Zos calculation using Hsig(Awave).
0.250 sz_alpha Factor for Wave dissipation surface tke flux .
100.000 crgban_cw Factor for Craig/Banner surface tke flux.

3.0000E-04 rdrg Linear bottom drag coefficient (m/s).
3.0000E-03 rdrg2 Quadratic bottom drag coefficient.
2.0000E-02 Zob Bottom roughness (m).
2.0000E-02 Zos Surface roughness (m).

2 Vtransform S-coordinate transformation equation.
4 Vstretching S-coordinate stretching function.
6.5000E+00 theta_s S-coordinate surface control parameter.
2.0000E+00 theta_b S-coordinate bottom control parameter.
100.000 Tcline S-coordinate surface/bottom layer width (m) used
in vertical coordinate stretching.

1025.000 rho0 Mean density (kg/m3) for Boussinesq approximation.
0.000 dstart Time-stamp assigned to model initialization (days).
20130101.00 time_ref Reference time for units attribute (yyyymmdd.dd)

3.6000E+02 Tnudg(01) Nudging/relaxation time scale (days)
for tracer 01: temp
3.6000E+02 Tnudg(02) Nudging/relaxation time scale (days)
for tracer 02: salt
3.0000E+01 Znudg Nudging/relaxation time scale (days)
for free-surface.
3.0000E+01 M2nudg Nudging/relaxation time scale (days)
for 2D momentum.
3.0000E+01 M3nudg Nudging/relaxation time scale (days)
for 3D momentum.
1.2000E+02 obcfac Factor between passive and active
open boundary conditions.

T VolCons(1) NLM western edge boundary volume conservation.
T VolCons(2) NLM southern edge boundary volume conservation.
T VolCons(3) NLM eastern edge boundary volume conservation.
T VolCons(4) NLM northern edge boundary volume conservation.

14.000 T0 Background potential temperature (C) constant.
35.000 S0 Background salinity (PSU) constant.
1.000 gamma2 Slipperiness variable: free-slip (1.0) or
no-slip (-1.0).

Thanks for your help!

[kate]

It looks like the standard pressure gradient error-driven currents. What are your Beckmann and Haney numbers? You didn't show that. Is your stratification a function of z only?

[Arvin_Yin]

It looks like the standard pressure gradient error-driven currents. What are your Beckmann and Haney numbers? You didn't show that. Is your stratification a function of z only?

Thank you, Kate!

My Beckmann and Haney numbers are
Maximum grid stiffness ratios:
rx0 = 3.661005E-01 (Beckmann and Haidvogel)
rx1 = 2.942017E+01 (Haney)

I have 54 vertical levels, and
! Set vertical, terrain-following coordinates transformation equation and
! stretching function (see below for details), [1:Ngrids].

Vtransform == 2 ! transformation equation
Vstretching == 4 ! stretching function

! Vertical S-coordinates parameters (see below for details), [1:Ngrids].

THETA_S == 6.5d0 ! surface stretching parameter
THETA_B == 2.0d0 ! bottom stretching parameter
TCLINE == 100.0d0 ! critical depth (m)

Stratification function is about Transformation Equations and Vertical Stretching Functions in Vertical S-coordinate page on ROMS Wiki? I never change it , so I think it is a function of z only.

And can you tell me more about the standard pressure gradient error? Thanks.

[kate]

Well, read the papers by Beckmann, Haney and Shchepetkin. Your rx numbers are similar to those I run with, but are considered quite large.

[Arvin_Yin]

Thanks.
I will read these papers. Then, what are the suitable range for rx numbers?

[jpringle]

When you are testing for the effects of pressure gradient error and topography, I would set all the boundaries to closed. I am worried about how the currents are reacting to both the northern and southern boundaries. I would not be surprised if the currents were substantially less if the boundaries were closed.

Essentially, open boundaries are hard. Pressure gradient errors are hard. Look at one at a time.

Jamie

[kate]

Then, what are the suitable range for rx numbers?

It's a trade-off. There is some value for which the model will blow up, so you want to smooth the bathymetry to get it below that. If you smooth too much, you'll be solving a significantly different problem and the results could be misleading in that way. For instance, Bernard Barnier once told me a story of a current system where the real world had a surface current with an undercurrent vertically stacked over each other. In the model with smoothing, the shelf break isn't as steep and there's room for one current to be beside the other, which changes the stability properties, etc.

So I run what I run, knowing there are some currents due to pressure gradient errors. One hopes those currents are significantly smaller than the wind-driven currents. It is just one of many sources of error, however, what with errors in forcing fields, boundary conditions, lack of resolution, etc.

P.S. Lack of resolution contributes to pressure gradient errors.

[rduran]

I would not be surprised if the currents were substantially less if the boundaries were closed

at the same time the signal can move around the domain when the boundaries are closed adding momentum each time it does, need to watch for that.


the rx parameters are not the whole story they are apriori upper bound estimates not sufficient conditions.

The following can help:

reduce number of sigma levels
decrease horizontal resolution,
smooth the bathymetry some more (specially where the error is biggest)
use a version of ROMS with a stiffened equation of state (I dont think ROMS rutgers has it and I dont know if getting it is a possibility)
combinations of the above

Edit: when I said "rx are apriori upper bound estimates not sufficient conditions" I meant the rx numbers recommended in the literature. It is not difficult to get considerable currents due to the pressure gradient with say r=0.2 (Beckman and Haidvogel r).

And regardless of whether your rx values are above or below the recommended values, it is not the end of the story, there are other things affecting the pressure gradient error not included in rx. The best way to know how much of an error you have is testing as described in the thread.

[Arvin_Yin]

Thank you for all the replies.
I will try one by one and report later.

[rduran]

Sorry I meant to say increase horizontal resolution, as Kate said.

Either decrease vertical resolution (e.g. Decreasing number of sigma levels) or increase horizontal resolution, or both.

You have to see what is required for your application, what error is acceptable and what error is not.