Unrealisitic thickness of mixed layer

[cyltby]

Dear all,


I am running a 3D ROMS of South China Sea with configuration of minimum horizontal grid size about 5 km and 30 vertical layers. Here HYCOM data products are applied as ROMS initial and boundary conditions and I choose MY2.5 for the vertical turbulent closure.

When I make a comparison of temperature profiles on two observation stations during Nov,2014, I find it inconceivable that ROMS mixed layer thickness is only about half of the observations while HYCOM performs much better.I also make a time series comparison of the temperature profile during the whole year of 2013. At that point, the HYCOM profile also match with the observation better than ROMS.

I believe our issue of the mixed layer depth should be some apparent error which may be due to misconfiguration.
I've attached all figures mentioned before and my head file&in file so that the situation would be easily understood.
Below are parts of the log file; Please let me know if there's any additional info I can provide to help explain this problem.

Since I am a novice ROMS users，any suggestion is highly appreciated.

### parts of the log file ###
Physical Parameters, Grid: 01
=============================

76199 ntimes Number of timesteps for 3-D equations.
415.000 dt Timestep size (s) for 3-D equations.
45 ndtfast Number of timesteps for 2-D equations between
each 3D timestep.
1 ERstr Starting ensemble/perturbation run number.
1 ERend Ending ensemble/perturbation run number.
12 nrrec Number of restart records to read from disk.
F LcycleRST Switch to recycle time-records in restart file.
6349 nRST Number of timesteps between the writing of data
into restart fields.
1 ninfo Number of timesteps between print of information
to standard output.
T ldefout Switch to create a new output NetCDF file(s).
6349 nHIS Number of timesteps between the writing fields
into history file.
1 ntsAVG Starting timestep for the accumulation of output
time-averaged data.
208 nAVG Number of timesteps between the writing of
time-averaged data into averages file.
1 ntsDIA Starting timestep for the accumulation of output
time-averaged diagnostics data.
6349 nDIA Number of timesteps between the writing of
time-averaged data into diagnostics file.
4.0000E+01 nl_tnu2(01) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 01: temp
4.0000E+01 nl_tnu2(02) NLM Horizontal, harmonic mixing coefficient
(m2/s) for tracer 02: salt
4.0000E+02 nl_visc2 NLM Horizontal, harmonic mixing coefficient
(m2/s) for momentum.
F LuvSponge Turning OFF sponge on horizontal momentum.
F LtracerSponge(01) Turning OFF sponge on tracer 01: temp
F LtracerSponge(02) Turning OFF sponge on tracer 02: salt
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.0000E-04 rdrg Linear bottom drag coefficient (m/s).
3.0000E-03 rdrg2 Quadratic bottom drag coefficient.
2.0000E-02 Zob Bottom roughness (m).
1.0000E+01 blk_ZQ Height (m) of surface air humidity measurement.
2.0000E+00 blk_ZT Height (m) of surface air temperature measurement.
1.0000E+01 blk_ZW Height (m) of surface winds measurement.
1 lmd_Jwt Jerlov water type.
2 Vtransform S-coordinate transformation equation.
4 Vstretching S-coordinate stretching function.
1.0000E+01 theta_s S-coordinate surface control parameter.
2.7000E+00 theta_b S-coordinate bottom control parameter.
50.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).
0.000 tide_start Reference time origin for tidal forcing (days).
20140101.00 time_ref Reference time for units attribute (yyyymmdd.dd)
1.0000E+00 Tnudg(01) Nudging/relaxation time scale (days)
for tracer 01: temp
1.0000E+00 Tnudg(02) Nudging/relaxation time scale (days)
for tracer 02: salt
1.0000E+00 Znudg Nudging/relaxation time scale (days)
for free-surface.
1.0000E+00 M2nudg Nudging/relaxation time scale (days)
for 2D momentum.
1.0000E+00 M3nudg Nudging/relaxation time scale (days)
for 3D momentum.
1.0000E-01 obcfac Factor between passive and active
open boundary conditions.
F VolCons(1) NLM western edge boundary volume conservation.
F VolCons(2) NLM southern edge boundary volume conservation.
F VolCons(3) NLM eastern edge boundary volume conservation.
F 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).
T LuvSrc Turning ON momentum point Sources/Sinks.
F LwSrc Turning OFF volume influx point Sources/Sinks.
F LtracerSrc(01) Turning OFF point Sources/Sinks on tracer 01: temp
T LtracerSrc(02) Turning ON point Sources/Sinks on tracer 02: salt
F LsshCLM Turning OFF processing of SSH climatology.
F Lm2CLM Turning OFF processing of 2D momentum climatology.
F Lm3CLM Turning OFF processing of 3D momentum climatology.
F LtracerCLM(01) Turning OFF processing of climatology tracer 01: temp
F LtracerCLM(02) Turning OFF processing of climatology tracer 02: salt
F LnudgeM2CLM Turning OFF nudging of 2D momentum climatology.
F LnudgeM3CLM Turning OFF nudging of 3D momentum climatology.
F LnudgeTCLM(01) Turning OFF nudging of climatology tracer 01: temp
F LnudgeTCLM(02) Turning OFF nudging of climatology tracer 02: salt
T Hout(idFsur) Write out free-surface.
T Hout(idu3dE) Write out 3D U-wastward component at RHO-points.
T Hout(idv3dN) Write out 3D V-northward component at RHO-points.
T Hout(idTvar) Write out tracer 01: temp
T Hout(idTvar) Write out tracer 02: salt

T Aout(idFsur) Write out averaged free-surface.
T Aout(idUvel) Write out averaged 3D U-momentum component.
T Aout(idVvel) Write out averaged 3D V-momentum component.
T Aout(idu3dE) Write out averaged 3D U-eastward at RHO-points.
T Aout(idv3dN) Write out averaged 3D V-northward at RHO-points.
T Aout(idWvel) Write out averaged W-momentum component.
T Aout(idTvar) Write out averaged tracer 01: temp
T Aout(idTvar) Write out averaged tracer 02: salt
T Aout(idDano) Write out averaged density anomaly.

[kate]

I was told probably ten years ago now that the MY2.5 scheme in ROMS is broken and should not be used. Instead, it is recommended to use either the GLS or the LMD schemes. Anyone know if the MY2.5 scheme was ever fixed? Or why it wasn't simply deleted?

[cyltby]

Hi Kate,

Thanks for your suggestion.Now I've used the LMD scheme instead and made another test.Since ROMS has blew up for many times,I have to reduce my time steps to make sure the model can go through.

Attachment is a new time series comparison of the temperature profile during the whole year of 2013（hot start from 2-1-2013）.Compared to the MY2.5 scheme which I posted before,the thickness of mixed layer remains unimproved.Maybe it is not the critical factor which vertical turbulence closure scheme that I select?

I am also wondering wheather I take my parameters for LMD scheme correctly or not. Any suggestion about parameterization for LMD scheme is appreciated.

[rduran]

I was told probably ten years ago now that the MY2.5 scheme in ROMS is broken and should not be used. Instead, it is recommended to use either the GLS or the LMD schemes. Anyone know if the MY2.5 scheme was ever fixed? Or why it wasn't simply deleted?

There may be researchers that are not fond of MY but a dedicated study concluded when comparing MY and KPP that

It would be useful if the differences found between the performance with the two parameterizations suggested one to be superior to the other. Indeed, in some aspects the additional physics included in the M-Y scheme, such as diffusion of TKE or buildup of TKE in a vertically confined boundary layer, seem intuitively more correct. However, a conclusive statement cannot be made without comparisons with real data and/or high-resolution numerical studies that resolve the fine-scale structure of these processes. Hopefully the points of difference noted in this study can help direct such future efforts.

http://marine.rutgers.edu/po/Papers/Durski_JGR04.pdf

Coastal modeling at OSU (including the author of that study) and many other well known researchers use MY.
-------

Re mixed layer depth, perhaps you need stronger wind or less incoming heat, are these values about right? Interpolating atmospheric forcing can be tricky.

[susonic]

Minor stuff: I guess you are using 4 times daily surface forcing and I don't think you need to use diurnal srflux in your cpp. Wind data is critical for your mixed layer depth calculation as Duran mentioned above.
BTW, is your SST reasonable compare with OBS?

My suggestion is
1. Use Akima 4th order advection on both horizontal and vertical.
if it's not working,
2. Use different vertical coordinate(say vtranform=2, vstretching=2)

Your log file will give us more information(whole log file).

-JH

[kate]

There may be researchers that are not fond of MY but a dedicated study concluded when comparing MY and KPP that

My comment was on the ROMS implementation of MY, not the algorithm itself. The MY2.5 implementation within the GLS routines was said to be more sound than the stand-alone routines.

As for getting the mixed layer depth right, it can be a function of many things, as susonic points out. Another thing to check is the Jerlov water type, which we now have as a function of depth, a sort of crude approximation of shallow waters being muddier.

[rduran]

My comment was on the ROMS implementation of MY, not the algorithm itself.

Thanks Kate I can see that now. It would be good to know if ROMS implementation has problems. At OSU ROMS is used quite a bit with MY, not sure what ROMS or MY for the matter.

Edit: In fact that paper uses ROMS and is over 11 years old

[cyltby]

Minor stuff: I guess you are using 4 times daily surface forcing and I don't think you need to use diurnal srflux in your cpp. Wind data is critical for your mixed layer depth calculation as Duran mentioned above.
BTW, is your SST reasonable compare with OBS?

My suggestion is
1. Use Akima 4th order advection on both horizontal and vertical.
if it's not working,
2. Use different vertical coordinate(say vtranform=2, vstretching=2)

Your log file will give us more information(whole log file).

-JH

Your comments on the above issues would be greatly appreciated！ It is necessary to double check my wind forcing file and SST. After that I will do the test cases as you suggested.
Sorry but, owing to my carelessness，the log file has been overwritten.

As for getting the mixed layer depth right, it can be a function of many things, as susonic points out. Another thing to check is the Jerlov water type, which we now have as a function of depth, a sort of crude approximation of shallow waters being muddier.

Thanks Kate for your reminder.What I am to do is to change Jerlov water type from type I to type IA,aiming to check out how it influences the mixed layer depth.

It helps me a lot.Thanks again!

[jzavala]

Hi Yong Chen,
We have a prototype model for the area that we have been fine-tuning and is now performing reasonably well.
The figure attached shows the temperature profile close to your station 30 and it compares well with your HYCOM plot.
Some specific aspects of our configuration relevant to your problem follow:
CPP options:
#define BULK_FLUXES
#define LONGWAVE_OUT
#define SOLAR_SOURCE
#define GLS_MIXING

In the input file we use the k-kl parameters for the GLS mixing, and water type = 1.
Our vertical discretization has 40 vertical levels, with Vtransform = 2, and Vstretching = 4.

For the surface fields required by BULK_FLUXES we use 3-hourly era-iterim analyses at the highest spatial resolution (1/8 deg).

If you need more specifics please feel free to send us a private email (jzavala@marine.rutgers.edu), we are always open for collaboration.

Javier

[johnluick]

I am responding to Kate's remark that the standalone MY mixing scheme may not be as reliable as the one in GLS. I did a simple test, which was two restarts of a well-tested coastal/shelf model, ran each for two years, saving five-day averages. One had the standalone, one had the GLS version. I looked at the MLD (as calculated by MLD_ecco.m by Katja Lorbacher and Dietmar Dommenget, which I think is as good as any for this purpose). There were seven CTD casts at a site in about 100 m depth). There was nothing in the results to say either mixing scheme was better than the other. The downwards evolution of the thermocline in autumn looked about the same in each, and the summer and winter minima and maxima were about the same. This isn't definitive but at least the standalone one is apparently not broken - this may be of some comfort to those who have been using it!

[rduran]

Thanks John. Your conclusion is similar to the paper I posted which also uses ROMS' implementations, that's good.

[rduran]

They found that KPP did better than MY2.5 for a coastal environment in this study, except during a downwelling episode -- see fig 9.

http://onlinelibrary.wiley.com/doi/10.1 ... 8/abstract

Circumstantial evidence is as good as it gets it seems, so the more of it the better.

[arango]

The MY2.5 scheme is coded correctly in ROMS and has several corrections that are not available in the original published formulation of the Mellor-Yamada 2.5. There are several papers in the literature about this. Also some people argue that the MY2.5 was not originally formulated correctly and the some of its approximations are incorrect... It is like we have two opposite views about this. Notice that the GLS scheme can be run as a MY2.5 scheme. See the parameters for the k-kl setup of the GLS. I had talked with Hans Burchard about this in the past. He argues that his formulation of the My2.5 in the GLS scheme solves several issues from the original Mellor and Yamada paper/code.

[chysun]

I am writing to respond Javier that the MLD in your temperature profile looks better than Yong Chen's. Yong Chen and I will do the test following your idea.Thanks!

[chysun]

Hi there,
I increase vertical level number from 30 to 40 and take k-kl parameters for GLS mixing instead folowing Javier's idea.Attachment shows that MLD changes little which indcates the problem relates to somthing else.I am gonna compare my forcing field with HYCOM forcing field if there is much difference.Anyway,the problem still remains unsolved.We sincerely look forward to any of your kindly reply!

[chysun]

Hi,
Good news on our issue! We check the tides at shore and found it unreasonably strong which might cause strong diffusion of thermocline at slope . The temp profile and time series make progress when we don't impose tide.

I guess we haven't set our tides correctly.We will double check our tide netcdf file.
some aspects of tidal setting
cpp file:
# !tidal setting
# define RADIATION_2D
# undef RAMP_TIDES
# define SSH_TIDES
# define UV_TIDES
# ifdef SSH_TIDES
# undef FSOBC_REDUCED
# define ADD_FSOBC
# endif
# if defined UV_TIDES || defined SSH_TIDES
# define ADD_M2OBC
# endif

Please correct me if there is anything wrong in the cpp file.


Best regards,
Haiyun