NEMO time-evolving sigma coordinates

[tiagoams]

Hi,

I am setting up parcels to run a 3D kernel using NEMO AMM7 CP6 mi-ax850. This is NEMO direct output as opposed to resampled fields that can be downloaded from CMEMS. The coordinates are time evolving sigma, with masked-out bottom cells to avoid excessive gradients.
I have uploaded a 3x3 subsection of the daily fields here: https://we.tl/t-HupV7snJa1

I am encountering several problems but I will start with the simplest: setting up the time evolving sigma coordinates, and what might be a silly mistake. I am following the tutorial and getting the error below.

Strangely the function is executed without errors until I replace
dimensions = {'U': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'depthw', 'time': 'time_counter'},
with
dimensions = {'U': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},

Your contributions would be greatly appreciated,
tiago

parcels.__version__

'2.4.0'

from parcels import FieldSet

data_path = './'

def prep_fieldset():#start,runtime, dt, outfname=None, outputdt=delta(hours=1)):#,kh_vertical=0):
    
    ufiles = sorted(glob(data_path+'20010101*grid_U.nc'))
    vfiles = sorted(glob(data_path+'20010101*grid_V.nc'))
    wfiles = sorted(glob(data_path+'20010101*grid_W.nc'))
    depthfiles = sorted(glob(data_path+'20010101*depthw.nc'))
    mesh_mask = data_path + 'coordinates_small.nc'

    filenames = {'U':      {'lon': mesh_mask, 'lat': mesh_mask, 'depth': depthfiles, 'data': ufiles},
                 'V':      {'lon': mesh_mask, 'lat': mesh_mask, 'depth': depthfiles, 'data': vfiles},
                 'W':      {'lon': mesh_mask, 'lat': mesh_mask, 'depth': depthfiles, 'data': wfiles},
                 'Kh_v3d': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': depthfiles, 'data': wfiles},
                }

    variables = {'U': 'vozocrtx',
                 'V': 'vomecrty',
                 'W': 'wo',
                 'Kh_v3d': 'avm',
                 'depthw': 'depthw',
                }

                                                        # should be deptht or depthw/depthu/depthv?
    dimensions = {'U':      {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'V':      {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'W':      {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'Kh_v3d': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'depthw': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                 }

    fieldset = FieldSet.from_nemo(filenames, variables, dimensions, allow_time_extrapolation=True)
    fieldset.U.set_depth_from_field(fieldset.depthw)
    fieldset.V.set_depth_from_field(fieldset.depthw)
    fieldset.W.set_depth_from_field(fieldset.depthw)
    fieldset.Kh_v3d.set_depth_from_field(fieldset.depthw)
    
    return fieldset

fieldset = prep_fieldset()

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~\AppData\Local\mambaforge\envs\py3_parcels\lib\site-packages\xarray\core\dataset.py:1319, in Dataset._construct_dataarray(self, name)
   1318 try:
-> 1319     variable = self._variables[name]
   1320 except KeyError:

KeyError: 'vozocrtx'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
Cell In [6], line 61
     57     fieldset.Kh_v3d.set_depth_from_field(fieldset.depthw)
     59     return fieldset
---> 61 fieldset = prep_fieldset()

Cell In [6], line 53, in prep_fieldset()
     45                                                     # should be deptht or depthw/depthu/depthv?
     46 dimensions = {'U':      {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
     47               'V':      {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
     48               'W':      {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
     49               'Kh_v3d': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
     50               'depthw': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
     51              }
---> 53 fieldset = FieldSet.from_nemo(filenames, variables, dimensions, allow_time_extrapolation=True)
     54 fieldset.U.set_depth_from_field(fieldset.depthw)
     55 fieldset.V.set_depth_from_field(fieldset.depthw)

File ~\AppData\Local\mambaforge\envs\py3_parcels\lib\site-packages\parcels\fieldset.py:494, in FieldSet.from_nemo(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, chunksize, **kwargs)
    492 if kwargs.pop('gridindexingtype', 'nemo') != 'nemo':
    493     raise ValueError("gridindexingtype must be 'nemo' in FieldSet.from_nemo(). Use FieldSet.from_c_grid_dataset otherwise")
--> 494 fieldset = cls.from_c_grid_dataset(filenames, variables, dimensions, mesh=mesh, indices=indices, time_periodic=time_periodic,
    495                                    allow_time_extrapolation=allow_time_extrapolation, tracer_interp_method=tracer_interp_method,
    496                                    chunksize=chunksize, gridindexingtype='nemo', **kwargs)
    497 if hasattr(fieldset, 'W'):
    498     fieldset.W.set_scaling_factor(-1.)

File ~\AppData\Local\mambaforge\envs\py3_parcels\lib\site-packages\parcels\fieldset.py:609, in FieldSet.from_c_grid_dataset(cls, filenames, variables, dimensions, indices, mesh, allow_time_extrapolation, time_periodic, tracer_interp_method, gridindexingtype, chunksize, **kwargs)
    606 if 'creation_log' not in kwargs.keys():
    607     kwargs['creation_log'] = 'from_c_grid_dataset'
--> 609 return cls.from_netcdf(filenames, variables, dimensions, mesh=mesh, indices=indices, time_periodic=time_periodic,
    610                        allow_time_extrapolation=allow_time_extrapolation, interp_method=interp_method,
    611                        chunksize=chunksize, gridindexingtype=gridindexingtype, **kwargs)

File ~\AppData\Local\mambaforge\envs\py3_parcels\lib\site-packages\parcels\fieldset.py:416, in FieldSet.from_netcdf(cls, filenames, variables, dimensions, indices, fieldtype, mesh, timestamps, allow_time_extrapolation, time_periodic, deferred_load, chunksize, **kwargs)
    414                     dFiles = fields[procvar].dataFiles
    415                     break
--> 416     fields[var] = Field.from_netcdf(paths, (var, name), dims, inds, grid=grid, mesh=mesh, timestamps=timestamps,
    417                                     allow_time_extrapolation=allow_time_extrapolation,
    418                                     time_periodic=time_periodic, deferred_load=deferred_load,
    419                                     fieldtype=fieldtype, chunksize=varchunksize, dataFiles=dFiles, **kwargs)
    421 u = fields.pop('U', None)
    422 v = fields.pop('V', None)

File ~\AppData\Local\mambaforge\envs\py3_parcels\lib\site-packages\parcels\field.py:372, in Field.from_netcdf(cls, filenames, variable, dimensions, indices, grid, mesh, timestamps, allow_time_extrapolation, time_periodic, deferred_load, **kwargs)
    370 filebuffer.name = filebuffer.parse_name(variable[1])
    371 if dimensions['depth'] == 'not_yet_set':
--> 372     depth = filebuffer.depth_dimensions
    373     kwargs['depth_field'] = 'not_yet_set'
    374 else:

File ~\AppData\Local\mambaforge\envs\py3_parcels\lib\site-packages\parcels\fieldfilebuffer.py:131, in NetcdfFileBuffer.depth_dimensions(self)
    128 @property
    129 def depth_dimensions(self):
    130     if 'depth' in self.dimensions:
--> 131         data = self.dataset[self.name]
    132         depthsize = data.shape[-3]
    133         self.data_full_zdim = depthsize

File ~\AppData\Local\mambaforge\envs\py3_parcels\lib\site-packages\xarray\core\dataset.py:1418, in Dataset.__getitem__(self, key)
   1416     return self.isel(**key)
   1417 if utils.hashable(key):
-> 1418     return self._construct_dataarray(key)
   1419 if utils.iterable_of_hashable(key):
   1420     return self._copy_listed(key)

File ~\AppData\Local\mambaforge\envs\py3_parcels\lib\site-packages\xarray\core\dataset.py:1321, in Dataset._construct_dataarray(self, name)
   1319     variable = self._variables[name]
   1320 except KeyError:
-> 1321     _, name, variable = _get_virtual_variable(self._variables, name, self.dims)
   1323 needed_dims = set(variable.dims)
   1325 coords: dict[Hashable, Variable] = {}

File ~\AppData\Local\mambaforge\envs\py3_parcels\lib\site-packages\xarray\core\dataset.py:176, in _get_virtual_variable(variables, key, dim_sizes)
    174 split_key = key.split(".", 1)
    175 if len(split_key) != 2:
--> 176     raise KeyError(key)
    178 ref_name, var_name = split_key
    179 ref_var = variables[ref_name]

KeyError: 'vozocrtx'

[erikvansebille]

Thanks for raising this issue, @tiagoams. I had a quick look at your code (thanks for providing the data for this minimal breaking example!) and think the difference is that in your case the depthw fields are in a separate file from the U, V , W etc fields. In the tutorial, all these fields are in the same file.

Is there a way in which you can merge your files so that the depthw coordinate is part of the U and V files too? Or that all variables are int he same file?

[tiagoams]

Thanks Erik, I have added the depth variable as "depthw4d" to one of the files and it now loads. Below is the code I used to prepare the fieldset and the modified forcing files are linked here

from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4, AdvectionRK4_3D, plotTrajectoriesFile,ErrorCode
#from parcels import ParcelsRandom, Field, DiffusionUniformKh, ErrorCode, Variable
from glob import glob
import numpy as np
from datetime import timedelta as delta
from datetime import datetime
from os import path
import xarray as xr
import matplotlib.pyplot as plt
import math


data_path = 'C:\\Users\\tams00\\Documents\\AMM7\\small\\'

def prep_fieldset():
    
    ufiles = sorted(glob(data_path+'20010101*_small_grid_U.nc'))
    vfiles = sorted(glob(data_path+'20010101*_small_grid_V.nc'))
    wfiles = sorted(glob(data_path+'20010101*_small_grid_W_depthw4d.nc'))
    mesh_mask = data_path + 'coordinates_small.nc'

    filenames = {'U':      {'lon': mesh_mask, 'lat': mesh_mask, 'depth': ufiles, 'data': ufiles},
                 'V':      {'lon': mesh_mask, 'lat': mesh_mask, 'depth': vfiles, 'data': vfiles},
                 'W':      {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': wfiles},
                 'Kh_v3d': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': wfiles},
                 'depthw4d': {'lon': mesh_mask, 'lat': mesh_mask, 'depth': wfiles, 'data': wfiles},
                }

    variables = {'U': 'vozocrtx',
                 'V': 'vomecrty',
                 'W': 'wo',
                 'Kh_v3d': 'avm',
                 'depthw4d': 'depthw4d',
                }

    dimensions = {'U':      {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'V':      {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'W':      {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'Kh_v3d': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                  'depthw4d': {'lon': 'glamf', 'lat': 'gphif', 'depth': 'not_yet_set', 'time': 'time_counter'},
                 }

    fieldset = FieldSet.from_nemo(filenames, variables, dimensions, allow_time_extrapolation=True)
    fieldset.U.set_depth_from_field(fieldset.depthw4d)
    fieldset.V.set_depth_from_field(fieldset.depthw4d)
    fieldset.W.set_depth_from_field(fieldset.depthw4d)
    fieldset.Kh_v3d.set_depth_from_field(fieldset.depthw4d)
    
    return fieldset

fieldset = prep_fieldset()

I started by doing a simple test by using AdvectionRK4 kernel, but I don't see the particle moving vertically with the free surface. Is this incorrect or am I missing something?

npart = 1
pset = ParticleSet.from_list(fieldset=fieldset, 
                             pclass=JITParticle,
                             lon=[-6.33]*npart,
                             lat=[49.93]*npart,
                             depth=np.linspace(3,3,npart)
                            )

kernels = pset.Kernel(AdvectionRK4)
pfile = pset.ParticleFile(data_path+"temp3.zarr", outputdt=delta(minutes=5))
pset.execute(kernels, runtime=delta(hours=2), dt=delta(minutes=5),  recovery={ErrorCode.ErrorThroughSurface: DeleteParticle}, output_file=pfile)

ds = xr.open_dataset(data_path+'temp3.zarr', engine='zarr')

fig, axs = plt.subplots(nrows=1, ncols=1, figsize=(6,2))

axs.plot(ds.time.T, ds.z.T,linewidth=1, marker ='.')
axs.invert_yaxis()
axs.set_ylabel('Depth(m)')
axs.set_xlabel('Date and time')
axs.set_title('Particle position')

On a second test, I tried to check where the top surface lies in the parcels referential. I released particles close to the surface and used a simple vertical diffusion kernel parameterised by diffusivity red from NEMO. The particles seem to be deleted at varying depths, which go tme confused. A few questions:

• should depth read from the files be the depth of water or distance below geoid (I might need to subtract the position of the free surface in the original files)?
• Is it correct to set the depth of variables in the U and V grid using the depthw varaible as I did above? i.e. does parcels correct for grid offset as in the horizontal grid?

Many thanks!

def VertDiffusionKh_simple(particle, fieldset, time):
    """Kernel for simple 3D diffusion where diffusivity (Kh) is provided by model (Kh_v3d).

    From DiffusionUniformKh:
    This kernel assumes diffusivity gradients are zero and is therefore more efficient.
    Since the perturbation due to diffusion is in this case isotropic independent, this
    kernel contains no advection and can be used in combination with a seperate
    advection kernel.

    The Wiener increment `dW` is normally distributed with zero
    mean and a standard deviation of sqrt(dt).
    """
    # Wiener increment with zero mean and std of sqrt(dt)
    dWz = ParcelsRandom.normalvariate(0, math.sqrt(math.fabs(particle.dt)))

    bz = math.sqrt(2 * fieldset.Kh_v3d[particle])    
    particle.depth +=  bz * dWz

    
# Passed to execute as recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle} 
def DeleteParticle(particle, fieldset, time):
    particle.delete()

# Release particles close to the surface
npart = 20
pset = ParticleSet.from_list(fieldset=fieldset, 
                             pclass=JITParticle,
                             lon=[-6.33]*npart,
                             lat=[49.93]*npart,
                             depth=np.linspace(0,3,npart)
                            )

kernels = pset.Kernel(VertDiffusionKh_simple)
pfile = pset.ParticleFile(data_path+"temp.zarr", outputdt=delta(seconds=5))
pset.execute(kernels, runtime=delta(minutes=5), dt=delta(seconds=5),  recovery={ErrorCode.ErrorThroughSurface: DeleteParticle}, output_file=pfile)

# Plot vertical position
ds = xr.open_dataset(data_path+'temp.zarr', engine='zarr')

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(10,4), sharey=True)

axs[0].plot(ds.time.T, ds.z.T)
axs[0].invert_yaxis()
axs[0].set_ylabel('Depth(m)')
axs[0].set_xlabel('Date and time')
axs[0].set_title('Particle position')

axs[1].plot(fieldset.depthw4d.data[:,0:4,1,1])
axs[1].set_xlabel('Time')
axs[1].set_title('Levels: depthw4d at 1,1')

[tiagoams]

The solution to my problem boils around answering the questions below, and doesn't require reading or running anything above :)

When using FieldSet.from_nemo:

• should depth in the input file be the depth below water surface or the distance below geoid (do I need to subtract the position of the free surface in the original files)?
• Is it correct to set the depth of variables in the U and V grid using the depth for the W grid (depthw variable)? i.e. does parcels correct for grid offset in the vertical as it does in the horizontal grid?

Many thanks

[erikvansebille]

should depth in the input file be the depth below water surface or the distance below geoid (do I need to subtract the position of the free surface in the original files)?

That I am not sure about. We normally don't do that in our simulations. I expect the differences between the two approaches to be very small, because it only matters for particles that are very close to the vertical cell boundaries. Since horizontal velocity interpolation is constant in the vertical within a cell for C-grids (see Delandmeter et al 2019), the only differences are when a particle is assumed to be in another cell. If the free surface variation is order centimeters, this happens in less than 0.1% for a 10m-think layer.

Is it correct to set the depth of variables in the U and V grid using the depth for the W grid (depthw variable)? i.e. does parcels correct for grid offset in the vertical as it does in the horizontal grid?

Yes, that is correct. Parcels needs to to know where the corner of the grid cell is, so it can do the right C-grid interpolation. Again, see Delandmeter et al 2019 for the reasons why

[tiagoams]

Many thanks Erik. I am concerned with getting the position of the surface and of the seabed correctly, so that particles moving verticaly due to vertical diffusivity can reflect at the surface or settle at the seabed (without the kernel crashing?). Please let me know if you have any suggestiong how to deal with this. I might have to use arbitrary buffer zones to avoid particles crossing the boundaries.

That I am not sure about. We normally don't do that in our simulations. I expect the differences between the two approaches to be very small, because it only matters for particles that are very close to the vertical cell boundaries. Since horizontal velocity interpolation is constant in the vertical within a cell for C-grids (see Delandmeter et al 2019), the only differences are when a particle is assumed to be in another cell. If the free surface variation is order centimeters, this happens in less than 0.1% for a 10m-think layer.

The free surface can move 10 meters in some coastal locations (e.g. Bristol Channel), or I am not understanding this correctly?

[erikvansebille]

The free surface can move 10 meters in some coastal locations (e.g. Bristol Channel), or I am not understanding this correctly?
Ah, I did not know that you were using Parcels in such an extreme condition. This clarifies your concerns.

I am concerned with getting the position of the surface and of the seabed correctly, so that particles moving verticaly due to vertical diffusivity can reflect at the surface or settle at the seabed (without the kernel crashing?). Please let me know if you have any suggestiong how to deal with this. I might have to use arbitrary buffer zones to avoid particles crossing the boundaries.

If you have problems with the RecoveryKernels, you could try using the new v3.0 branch of Parcels (#1426), which has a major reimplementation of error kernels; see https://github.com/OceanParcels/parcels/blob/v3.0/docs/examples/tutorial_kernelloop.ipynb.

[tiagoams]

My experiments show that the depth set in set_depth_from_field should be water depth and not position below mean sea level or geoid. I am getting a predictable behaviour from RecoveryKernels: ErrorThroughSurface when particle.depth<0 and ErrorOutOfBounds when particle depth is > depth of bottom layer (seabed layer), so I am able to catch those and nudge the particles in. Unfortunately ErrorOutOfBounds can mean that the particle goes out of bounds horizontally or through the seabed.

It is great to hear of the developments for version 3.0. Modifying particle_dlon/_ddepth instead of the position directly will be a great improvment for 3D diffusivity and animal movement. I will have to rethink my kernels when 3.0 is released.

[erikvansebille]

I will have to rethink my kernels when 3.0 is released.

You could already start using the 3.0 pre-release if you wanted, by installing the developer version and then checking out the v3.0 branch associated with #1426

[LauraBachmaier]

@tiagoams Is access to the NEMO direct output as opposed to resampled fields from Copernicus available publicly? @erikvansebille I am currently trying to set up a fieldset using NEMO data downloaded from Copernicus but it seems this data has been re-interpolated to an A-grid and I am unsure now whether the fieldset is being initiated correctly or not. Have you come across this issue? Any advice?

[JamiePringle]

Do not use the A-grid data. Because of the interpolation, there is velocity that goes through the land boundaries, leading to grounding of particles and other anomalies. Save yourself the several weeks my graduate student spent seeing just how bad the A-grid data is. Jamie

[LauraBachmaier]

@JamiePringle That's what I assumed and is exactly the issue I am encountering at the moment. Thank you for the advice!

[erikvansebille]

Thank you so much for the quick reply. So if I am using data through Copernicus and using fieldset.from_netcdf would this load the data in correctly (A-grid data loaded in using an A-grid method)? Or would this lead to interpolation problems?

No, you have to use the from_nemo method, that interpolates the C-grid correctly. Please refer to https://docs.oceanparcels.org/en/latest/examples/tutorial_nemo_3D.html

[tiagoams]

It is good to know that Mercator makes the original files of global GLORYS output available. For completeness, I am using AMM7 and AMM15 for the North West European Shelf developped by the Metoffice. To access raw outputs you either need to run them yourself, or request ftp access here. Another important difference of the raw output is the inclusion of vertical diffusivity in the output which is required to run a Visser (1997) type of kernel.

I am having several problems running with time-evolving sigma coordinates, but I will have to test with the latest version before posting details :)

[LauraBachmaier]

@tiagoams @erikvansebille Many thanks to both of you for your help on this! Great to know as well about the vertical diffusivity data inclusion. This is something I am also planning on including in the final version (though still some time away from this). I look forward to hearing more updates :)