Doubts regarding strong surface currents

[thiago.pdpaula]

Hi,

I'm running ROMS (svn 836) in a realistic application in southeast of Brazil and I've noticed what appears to be an overestimation of surface currents, possibly due to wind stress.

I'am using bulk fluxes and hourly fields from CFSR reanalyses as forcing and 1/12 degree resolution grid nested to Mercator Global, which provides initial conditions and daily boundary fields. These are my boundary conditions:

Code: Select all

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

 zeta                     1   Closed       Chapman Imp  Chapman Imp  Chapman Imp

 ubar                     1   Closed       Flather      Flather      Flather

 vbar                     1   Closed       Flather      Flather      Flather

 u                        1   Closed       Rad + Nud    Rad + Nud    Rad + Nud

 v                        1   Closed       Rad + Nud    Rad + Nud    Rad + Nud

 temp                     1   Closed       Rad + Nud    Rad + Nud    Rad + Nud

 salt                     1   Closed       Rad + Nud    Rad + Nud    Rad + Nud

 tke                      1   Closed       Gradient     Gradient     Gradient

And these are my CPP options:

Code: Select all

 Activated C-preprocessing Options:

 SSE                 ROMS SSE 112
 ANA_BSFLUX          Analytical kinematic bottom salinity flux.
 ANA_BTFLUX          Analytical kinematic bottom temperature flux.
 ASSUMED_SHAPE       Using assumed-shape arrays.
 BULK_FLUXES         Surface bulk fluxes parameterization.
 CURVGRID            Orthogonal curvilinear grid.
 DIFF_GRID           Horizontal diffusion coefficient scaled by grid size.
 DJ_GRADPS           Parabolic Splines density Jacobian (Shchepetkin, 2002).
 DOUBLE_PRECISION    Double precision arithmetic.
 EMINUSP             Compute Salt Flux using E-P.
 KANTHA_CLAYSON      Kantha and Clayson stability function formulation.
 LONGWAVE_OUT        Compute outgoing longwave radiation internally.
 MASKING             Land/Sea masking.
 MIX_GEO_TS          Mixing of tracers along geopotential surfaces.
 MIX_S_UV            Mixing of momentum along constant S-surfaces.
 MPI                 MPI distributed-memory configuration.
 MY25_MIXING         Mellor/Yamada Level-2.5 mixing closure.
 NONLINEAR           Nonlinear Model.
 NONLIN_EOS          Nonlinear Equation of State for seawater.
 N2S2_HORAVG         Horizontal smoothing of buoyancy and shear.
 POWER_LAW           Power-law shape time-averaging barotropic filter.
 PROFILE             Time profiling activated .
 K_GSCHEME           Third-order upstream advection of TKE fields.
 RADIATION_2D        Use tangential phase speed in radiation conditions.
 RI_SPLINES          Parabolic Spline Reconstruction for Richardson Number.
 !RST_SINGLE         Double precision fields in restart NetCDF file.
 SALINITY            Using salinity.
 SOLAR_SOURCE        Solar Radiation Source Term.
 SOLVE3D             Solving 3D Primitive Equations.
 SPLINES_VDIFF       Parabolic Spline Reconstruction for Vertical Diffusion.
 SPLINES_VVISC       Parabolic Spline Reconstruction for Vertical Viscosity.
 SPHERICAL           Spherical grid configuration.
 TS_C4HADVECTION     Fourth-order centered horizontal advection of tracers.
 TS_C4VADVECTION     Fourth-order centered vertical advection of tracers.
 TS_DIF2             Harmonic mixing of tracers.
 UV_ADV              Advection of momentum.
 UV_COR              Coriolis term.
 UV_U3HADVECTION     Third-order upstream horizontal advection of 3D momentum.
 UV_C4VADVECTION     Fourth-order centered vertical advection of momentum.
 UV_QDRAG            Quadratic bottom stress.
 UV_VIS2             Harmonic mixing of momentum.
 VAR_RHO_2D          Variable density barotropic mode.
 VISC_GRID           Horizontal viscosity coefficient scaled by grid size.

I set up a sponge with 15 grid points width in all boundaries with visc_factor and diff_factor respectively equal to 20 and 10. My current timestep is 180 seconds and my bathymetry was smoothed using Mathieu's Matlab package giving maximum rx0 and rx1 factors of 0.15 and 7.6, which I believe to be quite conservative.

The problem can be seen in the attached pictures, which shows a snapshot of surface currents and SSH from ROMS and from global MERCATOR model, to which I nested my grid, a few days after initialization. A snapshot of wind velocity is also shown.

I was expecting a geostrophic flow where currents vectors would follow more or less elevation gradients, as it happens on MERCATOR and in situ measurements show. However, currents on ROMS are stronger and presents a drift at the same regions where wind speed is also stronger, suggesting a possible overestimation of wind stress in my model.

At the beginning I thought this behavior would be restricted to the very near surface layer. However, the drift can be seen in deeper layers down to a ~ 50 m depth or more, specially on shallower regions.

The average flow seems to be overestimated as well, as shown in the pictures of averaged surface currents. ROMS currents are stronger than MERCATOR and other global models such as HYCOM NCODA, specifically along the slope in the domain of Brazil Current.

ROMS_average_surface_currents.png (32.14 KiB) Viewed 14560 times

MERCATOR_average_surface_currents.png (32.27 KiB) Viewed 14560 times

I made a few tests altering vertical mixing scheme (MY25, GLS k-epsilon, LMD), increasing viscosity (from 4.0 to 25.0 m2/s) and changing CFSR hourly data to ERA INTERIM 6 hourly fields, but I'm not satisfied yet. Near surface currents remains stronger than expected for the region and I wouldn't like to have mesoscale features vanished by high viscosity values.

A few doubts I have:
1. Are my high velocities and currents drift due to strong wind stress or have I misinterpreted my results? What else could be causing this behavior?
2. I heard that bulk fluxes may overestimate wind stress, specially on gust situations, so what parameters could I tweak in order to reduce wind stress and obtain more consistent values of surface currents?
3. Should I increase viscosity and/or reduce wind momentum components?
4. Should I switch any of my boundary conditions?
5. Are my configuration parameters set properly (see below)? Is there anything weird in my setup?

Other information about my run are summarized below.

Code: Select all

 Physical Parameters, Grid: 01
 =============================

     219360  ntimes            Number of timesteps for 3-D equations.
    180.000  dt                Timestep size (s) for 3-D equations.
         30  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.
          0  nrrec             Number of restart records to read from disk.
          T  LcycleRST         Switch to recycle time-records in restart file.
        480  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).
        240  nHIS              Number of timesteps between the writing fields
                                 into history file.
          0  nQCK              Number of timesteps between the writing fields
                                 into quicksave file.
 5.0000E+00  nl_tnu2(01)       NLM Horizontal, harmonic mixing coefficient
                                 (m2/s) for tracer 01: temp
 5.0000E+00  nl_tnu2(02)       NLM Horizontal, harmonic mixing coefficient
                                 (m2/s) for tracer 02: salt
 5.0000E+00  nl_visc2          NLM Horizontal, harmonic mixing coefficient
                                 (m2/s) for momentum.
          T  LuvSponge         Turning ON  sponge on horizontal momentum.
          T  LtracerSponge(01) Turning ON  sponge on tracer 01: temp
          T  LtracerSponge(02) Turning ON  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.
 1.0000E-02  Zob               Bottom roughness (m).
 2.0000E+00  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.
          3  lmd_Jwt           Jerlov water type.
          2  Vtransform        S-coordinate transformation equation.
          4  Vstretching       S-coordinate stretching function.
 5.0000E+00  theta_s           S-coordinate surface control parameter.
 3.0000E-01  theta_b           S-coordinate bottom  control parameter.
     40.000  Tcline            S-coordinate surface/bottom layer width (m) used
                                 in vertical coordinate stretching.
   1025.000  rho0              Mean density (kg/m3) for Boussinesq approximation.
  14153.000  dstart            Time-stamp assigned to model initialization (days).
19700101.00  time_ref          Reference time for units attribute (yyyymmdd.dd)
 3.0000E+01  Tnudg(01)         Nudging/relaxation time scale (days)
                                 for tracer 01: temp
 3.0000E+01  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.
 3.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).
          F  LuvSrc            Turning OFF 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
          F  LtracerSrc(02)    Turning OFF point Sources/Sinks on tracer 02: salt
          F  LsshCLM           Turning OFF processing of SSH climatology.
          F  Lm2CLM            Turning OFF processing of 2D momentum climatology.
          T  Lm3CLM            Turning ON  processing of 3D momentum climatology.
          T  LtracerCLM(01)    Turning ON  processing of climatology tracer 01: temp
          T  LtracerCLM(02)    Turning ON  processing of climatology tracer 02: salt
          F  LnudgeM2CLM       Turning OFF nudging of 2D momentum climatology.
          T  LnudgeM3CLM       Turning ON  nudging of 3D momentum climatology.
          T  LnudgeTCLM(01)    Turning ON  nudging of climatology tracer 01: temp
          T  LnudgeTCLM(02)    Turning ON  nudging of climatology tracer 02: salt

          T  Hout(idFsur)      Write out free-surface.
          T  Hout(idUbar)      Write out 2D U-momentum component.
          T  Hout(idVbar)      Write out 2D V-momentum component.
          T  Hout(idUvel)      Write out 3D U-momentum component.
          T  Hout(idVvel)      Write out 3D V-momentum component.
          T  Hout(idTvar)      Write out tracer 01: temp
          T  Hout(idTvar)      Write out tracer 02: salt
          T  Hout(idUsms)      Write out surface U-momentum stress.
          T  Hout(idVsms)      Write out surface V-momentum stress.
          T  Hout(idTsur)      Write out surface net heat flux.
          T  Hout(idTsur)      Write out surface net salt flux.
          T  Hout(idSrad)      Write out shortwave radiation flux.
          T  Hout(idLrad)      Write out longwave radiation flux.
          T  Hout(idLhea)      Write out latent heat flux.
          T  Hout(idShea)      Write out sensible heat flux.
          T  Hout(idEmPf)      Write out E-P flux.
          T  Hout(idevap)      Write out evaporation rate.
          T  Hout(idrain)      Write out rain rate.

...

 Metrics information for Grid 01:
 ===============================

 Minimum X-grid spacing, DXmin =  8.02032179E+00 km
 Maximum X-grid spacing, DXmax =  9.02554196E+00 km
 Minimum Y-grid spacing, DYmin =  8.02323505E+00 km
 Maximum Y-grid spacing, DYmax =  9.02110875E+00 km
 Minimum Z-grid spacing, DZmin =  1.11179822E-01 m
 Maximum Z-grid spacing, DZmax =  5.24095214E+02 m

 Minimum barotropic Courant Number =  7.89417705E-03
 Maximum barotropic Courant Number =  2.30942340E-01
 Maximum Coriolis   Courant Number =  1.30899694E-02

 Horizontal mixing scaled by grid size, GRDMAX =  9.02332508E+00 km

...

 Basin information for Grid 01:

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

 Initial basin volumes: TotVolume =  6.3595735250E+15 m3
                        MinVolume =  1.1017982474E+07 m3
                        MaxVolume =  4.0811095033E+10 m3
                          Max/Min =  3.7040442866E+03

Thanks in advance for any help.

[mainara]

Were you able to solve the problem?

[thiago.pdpaula]

I don't remember the exact steps I took, but the flag WIND_MINUS_CURRENT (see: https://www.myroms.org/projects/src/ticket/758) along with a fine tweak in viscosity (nl_visc2) led me to better results.

Besides comparing to global models, I also made some comparisons with in situ data (with appropriate time filtering) and the modeled velocities were in a reasonable range.

[mainara]

Thanks for your answer Thiago,

Would you answer me if you used 2 or 10 in blk_ZW and what was the value that you used for VISC2 or add_VISC2?

I'm studying the same region and I use:

blk_zw=10
VISC2=5000
add_VISC2=0

I don't remember the exact steps I took, but the flag WIND_MINUS_CURRENT (see: https://www.myroms.org/projects/src/ticket/758) along with a fine tweak in viscosity (nl_visc2) led me to better results.

Besides comparing to global models, I also made some comparisons with in situ data (with appropriate time filtering) and the modeled velocities were in a reasonable range.

[thiago.pdpaula]

Hi,

the BLK_ZW parameter refers to height of your surface wind. Mine is set to 10.0 m as this is the height of the wind momentum components I used to build my forcing files. I used this dataset:
"NCEP Climate Forecast System Reanalysis (CFSR) Selected Hourly Time-Series Products, January 1979 to December 2010" (https://rda.ucar.edu/datasets/ds093.1/).

Check the variable attributes of the datasets you used to build your forcing files to be sure about the correct value.

Regarding the horizontal viscosity, I think it really depends on your setup and grid resolution, thus, it is very unlike that my values will behave the same way for you. Anyway, in this application I used harmonic horizontal viscosity (UV_VIS2) with VISC2 == 20.0 m2/s.

You don't have to worry about ad_VISC2 if you're not running the Adjoint model.

[mainara]

Thanks for your quickly answer and help.

[hpd14thu]

Hello， I wonder that the magnitute of the wind speed you used.

I also met the problem that the currents were too strong. My wind speed is up tp 58m/s. I don't know that if the scheme in ROMS can adapt to the 58m/s-wind speed, or if there is something wrong in my cpps.

[mainara]

Are you using COAWST or ROMS?

Hello， I wonder that the magnitute of the wind speed you used.

I also met the problem that the currents were too strong. My wind speed is up tp 58m/s. I don't know that if the scheme in ROMS can adapt to the 58m/s-wind speed, or if there is something wrong in my cpps.

[hpd14thu]

Are you using COAWST or ROMS?

Hello， I wonder that the magnitute of the wind speed you used.

I also met the problem that the currents were too strong. My wind speed is up tp 58m/s. I don't know that if the scheme in ROMS can adapt to the 58m/s-wind speed, or if there is something wrong in my cpps.

Mainara, Thank you for your reply. I am using COAWST, but I just use the ROMS part in COAWST.I don't difine the WRF or SWAN

[jivica]

I suggest you first to read a few papers related to the wind stress formulations in extreme events like TC, many of them are John's papers using COAWST and many Hurricane cases etc. Not to just blindly run ROMS model and then expect to get result that is perfect.

For example:
Bryant & Akbar in JMSE 2016 doi:10.3390/jmse4030058 where you can see limits for drag formulations or
Hwang in JPO 2018 doi:10.1175/JPO-D-18-0107.1 for new formulation having more physical background and saturation of wind drag for TC at wind speeds > 30m/s

As you wrote, you are using COAWST, then it would make sense to use coupling and have better formulation coming from the wave model which is then fed back to WRF model providing forcing for ROMS model. This is the way to go, step by step.

Cheers,
Ivica

[mainara]

But, are you only using ROMS why use the COAWST ?
I agree with Ivica, if you do not use the functions of a coupled model, it does not make sense to use it.
Do you use the tidal forcing and river data?
What atmospheric and ocean forcing are you using?
On what machine do you perform your simulations?

I noticed that you did not use some CPPs, such as:
#define SPLINES_VDIFF
#define SPLINES_VVISC
#define VISC_GRID
# define DIFF_GRID
# define LMD_CONVEC
# define LMD_SKPP
# define LMD_NONLOCAL
# define LMD_SHAPIRO
# define ALBEDO_CURVE

I do not have much modeling expertise, but I can attachment to you my .in and .h in the study area.
If you want, mainarabg@gmail.com

[hpd14thu]

I suggest you first to read a few papers related to the wind stress formulations in extreme events like TC, many of them are John's papers using COAWST and many Hurricane cases etc. Not to just blindly run ROMS model and then expect to get result that is perfect.

For example:
Bryant & Akbar in JMSE 2016 doi:10.3390/jmse4030058 where you can see limits for drag formulations or
Hwang in JPO 2018 doi:10.1175/JPO-D-18-0107.1 for new formulation having more physical background and saturation of wind drag for TC at wind speeds > 30m/s

As you wrote, you are using COAWST, then it would make sense to use coupling and have better formulation coming from the wave model which is then fed back to WRF model providing forcing for ROMS model. This is the way to go, step by step.

Cheers,
Ivica

Thank you for your reply, Ivica. You really help me a lot. Wind stress is the issue point, and I will read more papers to figure out this.

Peida[/quote]

[hpd14thu]

But, are you only using ROMS why use the COAWST ?
I agree with Ivica, if you do not use the functions of a coupled model, it does not make sense to use it.
Do you use the tidal forcing and river data?
What atmospheric and ocean forcing are you using?
On what machine do you perform your simulations?

I noticed that you did not use some CPPs, such as:
#define SPLINES_VDIFF
#define SPLINES_VVISC
#define VISC_GRID
# define DIFF_GRID
# define LMD_CONVEC
# define LMD_SKPP
# define LMD_NONLOCAL
# define LMD_SHAPIRO
# define ALBEDO_CURVE

I do not have much modeling expertise, but I can attachment to you my .in and .h in the study area.
If you want, mainarabg@gmail.com

Thank you for your reply and I have found that the problem is the wind stress. After I modified the wind stress formula, I got better result.

Thank you!

Peida