Results and issues of barotropic tidal run

[BAN]

Hi all

I’m running ROMS in barotropic mode only forced by tides. The results are actually quite good for all validated constituents (surface and current) except for the current-ellipses of O1 in one area. In this area the modeled O1 current ellipses is up to five times too large (see Figure 1 below). The only other thing which seems suspicious is that the diurnal constituents seem to be somewhat out of sync with the forcing (see Figure 2 below). Beginning right inside the boundary, the modeled diurnal constituents (surface amplitude) are about one cm too large compared to the forcing and some boundary features are shifted somewhat. Similar behavior is not observed for the semi-diurnal constituents.

Here is what I have tried so far without getting rid of the problem:
Changing tunable parameters:
- Quadratic friction (rdrg) form 2*10^-3 to 3*10^-3
- Horizontal viscosity coefficient (visc2) from 200 to 20
Changing the forcing
- Run with surface forcing only
- Run with surface and current forcing
- Run with forcing from OSU tidal data inversion Atlantic model/domain
- Run with forcing from OSU tidal data inversion Global model/domain
Change the propagation dynamics
- Run with nonlinear advection terms
- Run without nonlinear advection terms

I’ve had some theories on what could be the cause, but now there are almost none left. Since this is my first ROMS application there is always the possibility that I’ve done some stupid mistake . . . see the first part of the log file below.

Question:
Is there any other parameter which could influence the results such as nudging/relaxation settings of Znudg, M2nudg, obcfac or others???

p.s. the model is first run for 10 days with any given settings before data are stored, so start-up effects should not be an issue.

Caption Figure 1:
Red lines are ROMS results (derived using t_tide). Black lines/dots are measured ADCP profiles.The picture in the bottom indicates depth in the model and the '+' is the location on the validation point. The abbreviations are maj (major axis), min (minor axis), inc (inclination), pha (phase) RCsp (residual current speed) and RCdir (direction of residual current). In each plot the crude average bias between model and data profile is given.

Caption Figure 2:
Uppermost plot is from the forcing O1 surface amplitude, the second plot is of the equivalent ROMS results. The bottom plot is the difference between the two. Colors indicate height in m. Black is land, black +’es are ADCP validation locations, dashed black lines are Lat and Long, and white lines/text give depth.

Beginning of log-file:
. . .
Physical Parameters, Grid: 01
=============================

Code: Select all

    1920000  ntimes          Number of timesteps for 3-D equations.
      1.800  dt              Timestep size (s) for 3-D equations.
          1  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.
       6000  nRST            Number of timesteps between the writing of data
                               into restart fields.
        200  ninfo           Number of timesteps between print of information
                               to standard output.
          T  ldefout         Switch to create a new output NetCDF file(s).
       1000  nHIS            Number of timesteps between the writing fields
                               into history file.
      24000  ndefHIS         Number of timesteps between creation of new
                               history files.
 2.0000E+01  nl_visc2        NLM Horizontal, harmonic mixing coefficient
                               (m2/s) for momentum.
 0.0000E+00  rdrg            Linear bottom drag coefficient (m/s).
 2.5000E-03  rdrg2           Quadratic bottom drag coefficient.
 0.0000E+00  Zob             Bottom roughness (m).
   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).
20100101.00  time_ref        Reference time for units attribute (yyyymmdd.dd)
 0.0000E+00  Znudg           Nudging/relaxation time scale (days)
                               for free-surface.
 0.0000E+00  M2nudg          Nudging/relaxation time scale (days)
                               for 2D momentum.
 2.0000E+01  obcfac          Factor between passive and active
                               open boundary conditions.
      1.000  gamma2          Slipperiness variable: free-slip (1.0) or 
                                                    no-slip (-1.0).
          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.

 Output/Input Files:

             Output Restart File:  faroegrid1_rst_nof.nc
        Prefix for History Files:  faroegrid1_his_nof
                 Input Grid File:  FaroeIslandsGrid1_grd.nc

 Tile partition information for Grid 01:  1403x0731  tiling: 004x001

     tile     Istr     Iend     Jstr     Jend     Npts

        0        1      351        1      731   256581
        1      352      702        1      731   256581
        2      703     1053        1      731   256581
        3     1054     1403        1      731   255850

 Tile minimum and maximum fractional grid coordinates:
   (interior points only)

     tile     Xmin     Xmax     Ymin     Ymax     grid

        0    -0.50  1404.50    -0.50   732.50  RHO-points
        1    -0.50  1404.50    -0.50   732.50  RHO-points
        2    -0.50  1404.50    -0.50   732.50  RHO-points
        3    -0.50  1404.50    -0.50   732.50  RHO-points

        0    -0.50  1404.50    -0.50   732.50    U-points
        1    -0.50  1404.50    -0.50   732.50    U-points
        2    -0.50  1404.50    -0.50   732.50    U-points
        3    -0.50  1404.50    -0.50   732.50    U-points

        0    -0.50  1404.50    -0.50   732.50    V-points
        1    -0.50  1404.50    -0.50   732.50    V-points
        2    -0.50  1404.50    -0.50   732.50    V-points
        3    -0.50  1404.50    -0.50   732.50    V-points

 Activated C-preprocessing Options:

 FAROEGRID1          2D-tidal model Faroe area grid 1 (grid: 1405x733, res.:500mx500m)
 ADD_FSOBC           Adding tidal elevation to proccesed OBC data.
 ANA_FSOBC           Analytical free-surface boundary conditions.
 ANA_INITIAL         Analytical initial conditions.
 ANA_M2OBC           Analytical 2D momentum boundary conditions.
 ANA_SMFLUX          Analytical kinematic surface momentum flux.
 ASSUMED_SHAPE       Using assumed-shape arrays.
 DOUBLE_PRECISION    Double precision arithmetic.
 EAST_FSCHAPMAN      Eastern edge, free-surface, Chapman condition.
 EAST_M2FLATHER      Eastern edge, 2D momentum, Flather condition.
 FSOBC_REDUCED       Using free-surface data in reduced physics condtions
 MASKING             Land/Sea masking.
 NONLINEAR           Nonlinear Model.
 NORTH_FSCHAPMAN     Northern edge, free-surface, Chapman condition.
 NORTH_M2FLATHER     Northern edge, 2D momentum, Flather condition.
 _OPENMP             OpenMP parallel shared-memory directives.
 PROFILE             Time profiling activated .
 RAMP_TIDES          Ramping tidal forcing for one day.
 !RST_SINGLE         Double precision fields in restart NetCDF file.
 SOUTH_FSCHAPMAN     Southern edge, free-surface, Chapman condition.
 SOUTH_M2FLATHER     Southern edge, 2D momentum, Flather condition.
 SSH_TIDES           Add tidal elevation to SSH climatology.
 UV_COR              Coriolis term.
 UV_QDRAG            Quadratic bottom stress.
 UV_VIS2             Harmonic mixing of momentum.
 WEST_FSCHAPMAN      Western edge, free-surface, Chapman condition.
 WEST_M2FLATHER      Western edge, 2D momentum, Flather condition.

 INITIAL: Configuring and initializing forward nonlinear model ...


 Minimum X-grid spacing, DXmin =  4.92730546E-01 km
 Maximum X-grid spacing, DXmax =  5.05848297E-01 km
 Minimum Y-grid spacing, DYmin =  4.92359303E-01 km
 Maximum Y-grid spacing, DYmax =  5.06168304E-01 km

 Minimum barotropic Courant Number =  3.36906469E-02
 Maximum barotropic Courant Number =  9.61345517E-01
 Maximum Coriolis   Courant Number =  2.37022567E-04

    GET_NGFLD   - tidal period
                   (Min =  4.30820453E+04 Max =  2.38068000E+06)
    GET_2DFLD   - tidal elevation amplitude
                   (Min =  0.00000000E+00 Max =  9.73707845E-01)
    GET_2DFLD   - tidal elevation phase angle
                   (Min =  1.14436078E-03 Max =  6.28303696E+00)

 Initial basin volumes: TotVolume =  1.2274133756E+14 m3
                        MinVolume =  1.1040633569E+06 m3
                        MaxVolume =  9.0342023445E+08 m3
                          Max/Min =  8.1826847053E+02

NL ROMS/TOMS: started time-stepping: (Grid: 01 TimeSteps: 00000001 - 01920000)

[/size]

[johnluick]

Ban

can you tell us some more about your adcp data, in particular - how many days went into the analysis behind the profiles?

John

[BAN]

Thank you for the interest John.

This particular ADCP deployment lasted just under a year, and the constituent information is supported by other deployments in the region (most are long time deployments). The data are from 150kHz RDI-ADCP and have undergone data quality procedures etc. The data have then been analyzed using the classic Forman Fortran code.

In other words I’m 99% sure the problem is in the model . . . or its forcing.

BAN

[johnluick]

BAN,

probably my arithmetic is wrong, but if I am correct, you only ran the model for 23 days or so, and since you don't use the first ten days, that would mean you only analysed only 13 days. If so I am surprised the results are so good!

John

[arango]

Yes, John is correct. Since you are using several tidal constituents, it will take a long simulation (perhaps, several months to a year) to resolve the harmonic decomposition of all these tidal constituents and the superposition of their frequencies. For example, around 30 days solution is needed to resolve the beat frequency of M2 and S2 for spring-neap tides. If you do the algebra, you will find out that 2*pi*(S2-M2)/2 ~ 29.57 days for just only 2 tidal components (M2=12.42h and S2=12.0h).

We needed like more tham 10 months of solution using ROMS option AVERAGES_DETIDE to resolve the harmonic least-square decomposition of 8 tidal constituents in the Philippine Archipelago. I really don't know what are you getting with just 13 days.

[BAN]

Hi John and arango, and thank you for your answers.

Hopefully there is some misunderstanding...
I believe that the model has been run for 40 days, and that I have excluded 10 of these as they most likely are contaminated by startup effects. This would leave the 30 days necessary in order to resolve neap spring tides.

Here is the arithmetic as I understand it:
Model timestep dt is 1.8 seconds
Number of timesteps ntimes is 1 920 000
This gives total time dt*ntimes = 3 456 000 seconds
Converted into days this is (3 456 000)/(60*60*24) = 40 days.
Current fields are written to history file every nHIS=1000 time steps, i.e. every 1800 seconds or 30 minutes.

History files have then been converted from model spatial snapshots into time series in each model point. These time series were then fed to t_tide telling it that the sampling was every 0.5 hours. t_tide was then used to derive resolved tidal constituents (with default Rayleigh factor=1) for every time series. In the t_tide results everything seems OK. S2, M2 etc. are resolved and are close to validation values. . . as you point out S2 and M2 would not be resolved if the time series was only 13 days and would at least be incomparable to reality if the reported sampling frequency was wrong.

It is true that I could have ‘cheated’ me to more frequencies than the length of the time series allows, by relaxing the Rayleigh restriction, but this is not done here . . . yet . . . according to Foreman and Henry (1989) this is OK as there is usually much less 'noise' in models than in reality. Turning off the advective terms in their model, they even resolve P1 and K1 (and K2 and S2) based only on a 30 day time series.

p.s.
If there is no clear mistake in the setup of the model, the reported O1-amplification problems could be a result of inappropriate boundary values. The ‘uppermost’ boundary of the model is from an area (Iceland-Faroe Ridge) which is known to contain other barotropic waves [Miller et al. (1996)]. These waves could perhaps (I have not got my head around the math yet) contaminate the daily OSU-constituents?

BAN

[rmueller]

How about bathymetry? Have you compared your ROMS model bathymetry to the actual values at the ADCP locations? Has smoothing (to meet Haney/rx0 criterion) altered the bathymetry in this location by an appreciable amount?

[BAN]

Thank you for the interest rmueller

Yes the depth at ADCP locations and the corresponding depth of the closest model point are mostly very similar. Take Figure 1 as an example the ADCP deployment estimated the depth to be 296 m, and the closest model point has 297 m depth.

Depth data have only been smoothed to meet the criteria of the eye, i.e. to smooth out the individual bottom tracks. I know that when running ROMS in 3d you must smooth bathymetry in order to alleviate ‘numerical currents’, but since this run is only in 2d, no smoothing has been done.

[rmueller]

Hi again,

Good point by you, re: the 2D bit. I will be interested to hear what you discover.

Best,
Rachael

[johnluick]

Ban,

did you resolve this issue? I've been hoping to hear you did. I think Hernan spelled it out clearly - and it wouldn't matter much whether you are talking about 13 days or 40 days. You need more days. Also, for a given number of days, you may get very different harmonic amplitudes depending on the number of constituents you analyze for. The least squares fit tries to apportion all the variance, so will apportion some to (say) O1 even if it really came from a neighboring frequency that was in the model but not in the analysis. You probably know all this. Personally, I would run the tidal model for a year, with all available TPXO constituents at the open boundary.

The reason it looks like an analysis problem is it is a diurnal component which is out of whack. Fewer oscillations, more time required to resolve.

John

[leonjld]

BAN, did you compare the ADCP currents directly with the tidal currents in the OSU tidal forcing? You can also try adcirc tidal database as forcing or compare with directly: http://www.unc.edu/ims/ccats/tides/tides.htm

if you have tidal elevation measurement it will be good to know first if the O1 elevation agrees or not.

[BAN]

Thank you for your remarks and questions John and leonjld

Re: John
- I’ll make longer time simulations later on, but I'm not sure this will make a big difference.

Re: leonjld
- As you can see below in Figure 3 the O1 overshoot is also present in OSU forcing data (green line is forcing, red is ROMS and black is data).

The solution must therefore be to find alternative forcing, at least for the daily constituents.

There is a difference (Fig 4 i.e. on the boundary) between forcing and ROMS, which I guess is caused by the fact that this run was forced by surface constituents only.

- I would love to use adcirc as forcing, but my domain lies in the Eastern North Atlantic and thus outside the coverage/range of adcirc.

- I have not validated against surface components yet, as only interest in the project is getting the currents right (project goal is tidal current energy estimates), but I’ll look into it.

QUESTION:
Is there any other tidal constituent database, with free access, that covers the East North Atlantic?

[ezaron]

Check out the FES2004 tidal model. It may be more accurate in this region.

-Ed

[BAN]

Thank you for the suggestion Ed.
I’ll look into it right away.

[BAN]

Hi again,
It is a long time since I last visited the forum. Here is a short update.

I tried out FES2004, and tried a lot of different model options, long time runs, fewer constituents etc. but the problem related to the unrealistically high O1 current speed north of the islands remained.

The procedure that reduced the problem most was to regenerate the grid (slightly coarser and with a non-tilted orientation) and rerun the whole thing. This gave O1 problems in the area which were compatible with the forcing model (OSU gave slightly better results). For some reasons, which I do not understand, the previous tilted grid amplified the O1 problem.

I also followed your advice and included surface measurements in the validation process.
The model actually performs good, both with respect to currents and surface, and I'm now working on a series of nested runs.

Some lessons learned along the way are:

Tuning of nl_visc2:
Has very little influence on the results, but lower values are slightly better and too low values make the model unstable.

Tuning of rdrg2:
Is quite important, as it has significant influence on the results. The specific value depends on model options. The best option in my specific case was to run with #undef UV_ADV and high friction (rdrg2 =0.0035). Almost as good was #define UV_ADV and lower friction (rdrg2 =0.0025).

Z vs. ZUV forcing
In these forums it is generally recommended to force the model by currents if possible, but in my specific case the results were better when only forced by surface. Another drawback of using current forcing (in large grids order of 1000 by 1500 points) is that is consumes over 40% of the computer time, (in smaller grids this is not a problem) in my specific case.

Use or non-use of #define FSOBC_REDUCED
The model does marginally better in runs forced by surface-only, when using FSOBC_REDUCED. This option does not cost much extra in computer time. It is nevertheless possible to get good results with #define SSH_TIDES only, but this is usually only possible if the time-step is reduced.

Use or non-use of #define UV_ADV
As mentioned above it is possible to get good results with this option turned off, but this makes it necessary to have higher friction in the model. The rationale seems to be that the energy that ought to go into into other tidal frequencies and the residual current (due to nonlinear propagation), can also be extracted by increased friction.