Time-varying depth in ROMS s-coordinate

[RomanIMMSP]

Dear Parcels community,

I need your help with the following issue.
I have ROMS model output to apply for particle tracking investigation. I studied the C-grid tutorial, time-varying depth tutorial, and discussion #1071, but problems still remain.
I have:
U-data with dimensions (time=7, depth=31, lat=560, lon=629)
V-data with dimensions (time=7, depth=31, lat=559, lon=630)
W-data with dimensions (time=7, depth=32, lat=560, lon=630)
I calculated changing in time depth
Z_w with dimensions (time=7, depth=32, lat=560, lon=630)
Z_rho with dimensions (time=7, depth=31, lat=560, lon=630)

Here is the code that I try to realize. I do not know if it assimilates U,V,W and depth data correctly. Probably someone with high experience in Parcels application can look through the code and give feedback. The code works but gives the message "ParticleSet is empty". Maybe the reason is in the depths. Initially, I set up depth 'depth': 's_rho' as s-values (from -1 to 0). Later I use the function fieldset.U.set_depth_from_field(fieldset.zw) to set up correct depth. Maybe it doesn't work correctly. I don't know how to check this. When I put 'depth': 'not_yet_set' it gives another error (KeyError: 'z_rho').

from parcels import AdvectionRK4_3D
from parcels import FieldSet
from parcels import JITParticle
from parcels import ParticleSet
from parcels import Variable
from parcels import plotTrajectoriesFile
from parcels import ErrorCode

import numpy as np
import math
import xarray as xr
from datetime import timedelta as delta


filenames = {'U': {'lon': 'Input_data/waom10extend_grd.nc', 'lat': 'Input_data/waom10extend_grd.nc', 'depth':'Input_data/ocean_avg_0001.nc', 'data': 'Input_data/ocean_avg_0*.nc'}, 
             'V': {'lon': 'Input_data/waom10extend_grd.nc', 'lat': 'Input_data/waom10extend_grd.nc', 'depth':'Input_data/ocean_avg_0001.nc', 'data': 'Input_data/ocean_avg_0*.nc'},
             'W': {'lon': 'Input_data/waom10extend_grd.nc', 'lat': 'Input_data/waom10extend_grd.nc', 'depth':'Input_data/ocean_avg_0001.nc', 'data': 'Input_data/ocean_avg_0*.nc'},
             'zw': {'lon': 'Input_data/waom10extend_grd.nc', 'lat': 'Input_data/waom10extend_grd.nc', 'depth':'Input_data/ocean_avg_0001.nc', 'data': 'Input_data/depth_rho_0*.nc'}}
variables = {'U': 'u', 'V': 'v', 'W': 'w', 'zw': 'z_rho'}
indices = {'U': {'lon': range(0, 629), 'lat': range(0, 559)},
           'V': {'lon': range(0, 629), 'lat': range(0, 559)}}
dimensions = {'U': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 's_rho', 'time': 'ocean_time'},
              'V': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 's_rho', 'time': 'ocean_time'},
              'W': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 's_rho', 'time': 'ocean_time'},
              'zw': {'lon': 'lon_rho', 'lat': 'lat_rho', 'depth': 's_rho', 'time': 'ocean_time'}}

fieldset = FieldSet.from_c_grid_dataset(filenames, variables, dimensions, gridindexingtype='mitgcm')

fieldset.U.set_depth_from_field(fieldset.zw)
fieldset.V.set_depth_from_field(fieldset.zw)
fieldset.W.set_depth_from_field(fieldset.zw)


pset = ParticleSet.from_line(fieldset=fieldset, pclass=JITParticle, size=31, start=(-25, -70), finish=(-25, -73), depth=-10)

def DeleteParticle(particle, fieldset, time):
    particle.delete()

output_file = pset.ParticleFile(name="output.nc", outputdt=delta(hours=6))

pset.execute(AdvectionRK4_3D, runtime=delta(days=20), dt=delta(minutes=5), output_file=output_file, recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})
pset.show()

output_file.export()
plotTrajectoriesFile('output.nc')

[erikvansebille]

Hi @RomanIMMSP, thanks for your question. I never work with time-varying depth coordinates myself, so have little experience. Perhaps @cjongedijk can help, she works with this type of data more often?

However the KeyError you mentioned suggests that z_rho simply doesn't exist in your netcdf files. May it be that you have to use s_rho instead of z_rho in your variables dictionary, and then 'depth': 'not_yet_set' in the dimensions dictionary, like in the tutorial?

[RomanIMMSP]

Hi erikvansebille.
Thanks for your prompt reply and suggestion. s_rho does not change in time and shows only vertical s-levels. It is 1D variable. While z_rho contains real depths for every point of the computation grid. It changes in time and has 4 dimensions, like U,V,W. Probably I should set up 'depth': 's_rho' as dimension for z_rho and 'depth': 'not_yet_set' for U,V,W. The part of code is below. When I made such changes, a new error appears (on the screen below). As I understand, there is a conflict between 32 vertical layers in W variable and 31 vertical layers in other variables. How can I solve this? If @cjongedijk can give her feedback, I will be very appreciated.

filenames = {'U': {'lon': 'Input_data/waom10extend_grd.nc', 'lat': 'Input_data/waom10extend_grd.nc', 'depth':'Input_data/ocean_avg_0001.nc', 'data': 'Input_data/ocean_avg_0*.nc'}, 
             'V': {'lon': 'Input_data/waom10extend_grd.nc', 'lat': 'Input_data/waom10extend_grd.nc', 'depth':'Input_data/ocean_avg_0001.nc', 'data': 'Input_data/ocean_avg_0*.nc'},
             'W': {'lon': 'Input_data/waom10extend_grd.nc', 'lat': 'Input_data/waom10extend_grd.nc', 'depth':'Input_data/ocean_avg_0001.nc', 'data': 'Input_data/ocean_avg_0*.nc'},
             'z_rho': {'lon': 'Input_data/waom10extend_grd.nc', 'lat': 'Input_data/waom10extend_grd.nc', 'depth':'Input_data/ocean_avg_0001.nc', 'data': 'Input_data/depth_rho_0*.nc'}}
variables = {'U': 'u', 'V': 'v', 'W': 'w', 'z_rho': 'z_rho'}
indices = {'U': {'lon': range(0, 629), 'lat': range(0, 559)},
           'V': {'lon': range(0, 629), 'lat': range(0, 559)}}
dimensions = {'U': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'not_yet_set', 'time': 'ocean_time'},
              'V': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'not_yet_set', 'time': 'ocean_time'},
              'W': {'lon': 'lon_psi', 'lat': 'lat_psi', 'depth': 'not_yet_set', 'time': 'ocean_time'},
              'z_rho': {'lon': 'lon_rho', 'lat': 'lat_rho', 'depth': 's_rho', 'time': 'ocean_time'}}

fieldset = FieldSet.from_c_grid_dataset(filenames, variables, dimensions, gridindexingtype='mitgcm')

fieldset.U.set_depth_from_field(fieldset.z_rho)
fieldset.V.set_depth_from_field(fieldset.z_rho)
fieldset.W.set_depth_from_field(fieldset.z_rho)

pset = ParticleSet.from_line(fieldset=fieldset, pclass=JITParticle, size=31, start=(-25, -70), finish=(-25, -73), depth=-10)

[erikvansebille]

Yes, this suggests that indeed the not_yet_set now works.

You could try to use indices={'W': {'depth': range(0, 30)}} (not tested) to force Parcels to read only the first 31 layers of the W variable? Although this may have unexpected consequences for the boundary conditions, so check carefully that the results make sense!

[RomanIMMSP]

I made this including, but nothing changed. I have the same error.

[erikvansebille]

Oh strange. And if you add indices={'depth': range(0, 10)}to the FieldSet creation?

[RomanIMMSP]

Hi Erik,

Do you mean in such way:
fieldset = FieldSet.from_c_grid_dataset(filenames, variables, dimensions, gridindexingtype='mitgcm', indices={'depth': range(0,30)})
In this case it works, but gives the message 'ParticleSet is empty'.

When I check the output file, I see 31 particles but only 2 time moments. Example:
traj obs time
"1 of 31" "1 of 2" 0.0
"1 of 31" "2 of 2" 0.0
"2 of 31" "1 of 2" 0.0
"2 of 31" "2 of 2" 0.0
"3 of 31" "1 of 2" 0.0
"3 of 31" "2 of 2" 900.0
"4 of 31" "1 of 2" 0.0
"4 of 31" "2 of 2" 900.0
Maybe the problem is in time interpretation. I have the following message, which, however, could be not important because in 2D case I had the same message and received trajectories for all particles (without message 'ParticleSet is empty').

Thanks for your support and have a good weekend.

[RomanIMMSP]

Hi @erikvansebille,

I made several tests during the weekend, and I am not sure that 'depth': 'not_yet_set' and set_depth_from_field work well in my code. Could you recommend how can I test/debug/check them?
In the current version of the code, I received the OutOfBounds Error immediately after the start of calculations. This error appears even when I apply fieldset.W.set_scaling_factor(0.) and change the initial depth of released particles.
I am not sure that it is important for Parcels, but in the ROMS model, the vertical grid indexing begins from the bottom and ends at the surface. Depth has negative values. Just as an example: depth(0)=-4500, depth(30)=-4.

Best wishes,
Roman

[erikvansebille]

Ah, it could very well be that this switched depth ordering creates challenges for the Parcels code. I have never seen such a depth structure before, so have no experience with how to deal with it in Parcels. Perhaps others here on Github have?

[RomanIMMSP]

Yes, you are right. The problem is in inverse indexing in the vertical direction. Since I am calculating z_rho separately, I change the indexing during the calculation using such a function depth11=depth1.reindex(s_rho=depth1.s_rho[::-1]). Now the example of depth in the z_rho variable looks like depth(0)=-4, depth(30)=-4500.
After such changes, the model starts to work without errors. But the results are incorrect because I also need to reindex U,V,W variables through the z-axis. Is it possible to do this within fieldset? Or maybe you can recommend some other method?

[RomanIMMSP]

Hi @erikvansebille,

I reindexed all needed variables outside Parcels as preprocessing to avoid inverse indexing in the vertical direction. Now the model works, hopefully without mistakes. I have several questions about model features.

1. In discussion strange behaviours with AdvetionRK4_3D in s-coord #765 you suggested using alternative Advection_3D kernels to describe the vertical movement of the particles in s-coordinates correctly. Is it applicable in my case?
2. Vertical positions of (u,v) and w are different (z_rho and z_w) on the C-grid. Which one should I use in the FieldSet definition? As I understand from the 'Nemo 3D tutorial' I should use z_w for all three velocity components. Is it true?
3. I have the situation when particles cross the surface and delete due to OutOfBoundError. According to my understanding, particles should be reflected from the surface back to the water. Does Parcels have such a feature or how to realize this?

Best wishes,
Roman

[erikvansebille]

Dear @RomanIMMSP; I'm sorry but I can't answer these questions. I am not working with ROMS myself, and these are choices that depend on your model setup and research question. You will have to answer your questions yourself

[daanreijnders]

About your third point: I've dealt with something similar as well. You can write a kernel that checks whether a particle is above a certain height and, if it crosses that depth, places is back:

def keep_below_depth(particle, fieldset, time):
    if particle.depth > -0.5:
        particle.depth = -0.5 # or whatever depth you use

Then you can simply chain this kernel after the advective step. For kernel chaining, see for instance this tutorial: https://nbviewer.org/github/OceanParcels/parcels/blob/master/parcels/examples/tutorial_sampling.ipynb?flush_cache=true

[RomanIMMSP]

Dear @daanreijnders, the problem is that position of the surface z_w(t,0,j,i) is different at every point of the modeling area and it changes in time. So, it should look like

if particle.depth > z_w(t,0,j,i):
        particle.depth = z_w(t,0,j,i)-0.1 # or something like that

I think that another solution could be in replacing all positive vertical velocities at the surface with the negative values. Then velocities near the surface will push particles back to the water.