Question: how does ROMS work for 3D without zeta included?

[kthyng]

I've been using the ROMS reader a lot and some of it has remained opaque to me. One question is: how is depth handled given that zeta (sea surface height) is not used and zeta is required to calculated the depth of the vertical layers in time for ROMS? Is it assuming quiescent surface for calculating the vertical layers?

[knutfrode]

I have in fact not used that reader much myself, but made the original version based on some local ROMS files which did not even contain zeta, as bathymetry is generally quite steep off the Norwegian coastline, so that zeta << H.

The sigma-to-z conversion is made using a computationally efficient method from the package Roppy of Bjørn Aadlandsvik (https://opendrift.github.io/_modules/opendrift/readers/roppy/depth.html#sdepth) where zeta is not taken into account. If zeta exists in the file, it will be returned as a regular variable (along with currents, temperature etc), that the "module developer" can access and use e.g. within the update method.

However, I am not sure whether this is strictly correct, or if zeta should have been taken into consideration in the conversion from sigma to z, so that z would be 0 at the very surface (zeta)? If so, I guess zeta should again be subtracted within the reader to be consistent with the general convention in OpenDrift that z refers to mean sea level. In other words, I guess the main question is whether zeta should be included in the stretching of vertical coordinates, or if it can be added afterwards? I would assume any difference would be negligible for our local applications where bathymetry is steep off the coast (zeta<<H), but it might perhaps be important for your applications that even include wetting and drying? @johannesro might have a comment to this.

(As perhaps already meantioned: it is a future goal that the generic netCDF-reader should also support other vertical coordinates than z, so that a dedicated ROMS reader should not be needed. Though that might require ROMS output to be strictly adhering to CF-convention, which it was not at the time when the ROMS reader was orinally made.)

[kthyng]

I have to stare a little harder at the roppy code to be sure but it looks the same as e.g. what I have in xroms but without zeta:

https://github.com/xoceanmodel/xroms/blob/c42b795e6e08cd77e038e2845c627d858c7705df/xroms/xroms.py#L202-L205

With this calculation, z is 0 at mean sea level (not the very surface except when the surface is instantaneously at mean sea level) and zeta varies positively and negatively from there. In the special case of surface-only simulations or when not available, zeta wouldn't be necessary.

whether zeta should be included in the stretching of vertical coordinates, or if it can be added afterwards? zeta is part of the definition of the vertical coords (though of course like you say some places it just doesn't matter).

It does sound like there is a difference in convention between what you've been doing in opendrift and what is typical in ROMS — I had wondered about this in the plotting I was just starting to do.

I will definitely need zeta included for the tidal flat areas I am working in.

It's hard for me to tell if I should put effort into the ROMS reader or getting my ROMS output to work in the generic reader. I am not sure if brand new ROMS output is CF-compliant or not yet, but I also still work with older output.

[knutfrode]

Making the generic reader support ROMS native output would be a long term goal, but I believe that a general adaptation of OpenDrift for better support of water and particles above z=0 should hopefully be fairly straightforward to do on short time, although it is a quite fundamental change.

It seems we then agree that z=0 should still mean the mean sea level, and that water and particles can then have z-values between seafloor_depth_below_mean_sea_level (constant in time) and zeta (variable in time).

I can think of a few necessary adaptations:

• Surface stokes drift (e.g. from wave model or parameterised) is presently applied fully where z=0, and declining with depth. But this should then be applied to z=zeta, and declining from there. The same goes for the wind_drift_factor.
• Presently turbulent mixing is reflecting particles at z=0, but this must now also happen at z=zeta.
• Thus zeta must now be always available as an environment variable, but will default to 0 if not provided e.g. from a ROMS file.
• The sigma-to-z convertion in ROMS reader might need to be updated, as I guess sigma=0 should correspond to z=zeta, and not z=0 as presently(?)
• probably some more changes here and there, but hopefully no major re-writing necessary.

I have some motivation to contribute to this, as we might also for the first time need this feature in a potential near future project.

[kthyng]

Oof, I was hoping the changes would only need to occur in the ROMS reader. I am happy to handle the changes in the ROMS reader but I am worried that my lack of familiarity in the other parts will leaves the changes piecemeal. Also, I'm planning to use multiple models (OpenOil, OceanDrift, LarvalFish, and Leeway) — will changes need to occur it these and others?

I could start on the ROMS reader tomorrow morning when I am fresh.

[knutfrode]

I believe changes would mainly be needed in OceanDrift module, which handles most of the generic 3D-things, as well as some of the methods from physics_methods, e.g. for Stokes drift. OpenOil, LarvalFish and PlastDrift etc are subclasses of OceanDrift, and little or no additional updates would be needed in those. Leeway is a 2D module, and should not be affected at all.

It could be a good share of work if you could try to improve the ROMS reader, as you know ROMS and sigma-coordinates better than myself. And then I could take care of the updates needed elesewhere. I assume you would like to see this implemented as soon as possible?

[kthyng]

I assume you would like to see this implemented as soon as possible?

Ideally, yes. If you're able, when might you be able to work on it? I've been using opendrift in the meantime just not realizing this issue is making it less accurate for me so in that sense I can keep using it for awhile like this but would nice to make changes and move on if possible.

[knutfrode]

There are some other urgent tasks as well, but I can try to get an overview tomorrow. It could probably be realistic to have this implemented within next week.

[kthyng]

Thus zeta must now be always available as an environment variable, but will default to 0 if not provided e.g. from a ROMS file

Does this mean that if zeta is not found in the ROMS file it should be assumed to be a scalar 0 (like a fallback value)? I'm trying to decide for how to write the code in depth.py.

[knutfrode]

Does this mean that if zeta is not found in the ROMS file it should be assumed to be a scalar 0 (like a fallback value)? I'm trying to decide for how to write the code in depth.py.

Yes, each reader does only know its own data, and the request it received: which variables, at given coordinates, depth and time. But OpenDrift will only request the variables that readers can actually provide. So if the ROMS dataset does not contain zeta, it will also never be requested to return zeta.

It is then up the the given simulation object and its configuration whether fallback zeta should be 0 or something else.

[kthyng]

@knutfrode I am opening up a PR so you can see what I've done and get feedback on next steps. You could work on the same PR if you want, also.

Other comments and questions:

1. Should the default in depths.py in the sdepth function for zeta be 0 or have the default for zeta be set elsewhere?
2. In depths.py in the sdepth function I allowed for zeta to be an array or a scalar — is this ok?
3. Are there other places I should add zeta in the ROMS reader from your understanding or so you would have access in other functions?

4. Thus zeta must now be always available as an environment variable, but will default to 0 if not provided e.g. from a ROMS file.

   Should I do something to make this available?

[knutfrode]

It is not clear to me how or if Roppy handles zeta different from 0.

Maybe the creator of Roppy is available to comment? @bjornaa: do you remember if or how Roppy deals with zeta different from 0? I do not see zeta explicitly mentioned anywhere in the Roppy code: https://github.com/bjornaa/roppy/blob/master/roppy/depth.py Does that mean that sigma=0 is always mapped to z=0? And if so, is this an approximation, or is it formally correct, and so that zeta should be added to z afterwards, so that z=zeta>0 actually corresponds to sigma<0?

Note that in the ROMS reader, the convenience method multi_zslize is used to interpolate the 3D-variables to some hardcoded z-levels for each pixel: https://github.com/bjornaa/roppy/blob/master/roppy/depth.py#L210 https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/reader_ROMS_native.py#L129

@kthyng The above interpolation to fixed z-levels might seem a bit awkward, but it is a consequence of OpenDrift being generic wrt any ocean or atmospheric models, and OpenDrift internally always use z as vertical coordinate, so that any reader is responsible to interpolate its own data to z-levels. The above fixed levels may though be overridden by the user. For a trajectory model written specicially for ROMS, these things would be simpler and more accurate, at the cost of the flexibility to easily combine forcing from (m)any types of models.

[kthyng]

I'm not sure I understand your question. If it is in response to my question about the default, I just meant I wasn't sure where to put the default of zeta of 0 (as in whether to have it in the sdepth function itself or at a higher level).

zeta isn't in the depth code you link in roppy to since it's assuming that it's not important relative to the other terms, as you've mentioned before. You can see how it should be in there in the code I implemented or in the links in the code I put in (from the CF definitions).

Does that mean that sigma=0 is always mapped to z=0? And if so, is this an approximation

This is an approximation.

[kthyng]

@knutfrode So if z is at depth it incorporates zeta

https://github.com/kthyng/opendrift/blob/54201a6f8805d565f72f23a024deaccea6cf7d62/opendrift/readers/reader_ROMS_native.py#L425-L434

However if it is at the surface it is not incorporated yet. Perhaps I could have variables['z'] be equal to zeta here?

https://github.com/kthyng/opendrift/blob/54201a6f8805d565f72f23a024deaccea6cf7d62/opendrift/readers/reader_ROMS_native.py#L414-L416

or do you have a different plan?

[knutfrode]

@knutfrode So if z is at depth it incorporates zeta https://github.com/kthyng/opendrift/blob/54201a6f8805d565f72f23a024deaccea6cf7d62/opendrift/readers/reader_ROMS_native.py#L425-L434

Ok, this looks good.

However if it is at the surface it is not incorporated yet. Perhaps I could have variables['z'] be equal to zeta here? https://github.com/kthyng/opendrift/blob/54201a6f8805d565f72f23a024deaccea6cf7d62/opendrift/readers/reader_ROMS_native.py#L414-L416 or do you have a different plan?

Until now, z=0 was always the ocean surface, so that z.min()==0 in this context means that all particles are at the surface (since z is here the depths of the actual particles in the request from outside (yes, I should have used more explicit variable names here) variables is the dictionary to return, and z is its vertical coordinate, which in this case was just 0, as all particles in this case were at surface, and just the surface layer from ROMS was returned.

But in the (near) future, I believe we must allow for particles in the ocean with positive z (previously assumed to be in air), whenever zeta is positive? Thus I guess z.min()==0 should be changed to z.min()>=0. However: variables is the dictionary to return with data, and z is its vertical coordinate(s), which have to be from this hardcoded array: https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/reader_ROMS_native.py#L129 This we cannot simply set variables["x"] to zeta, as zeta is varying over the domain, and can take "any" value.

I believe the reason for this and other things will be very difficult to understand without knowing the details of the relationship between OpenDrift and "Readers", thus an explanation of that is needed:

The ROMS reader, as any other readers, is in fact independent of the OpenDrift trajectory modules, and we have actually considered extracting the readers and their functionality (collected in class Environment) out of the OpenDrift package. The readers have some clearly defined tasks ("API"):

• Upon initialisation (__init__) the reader should open the files/URLs/... and detect the coverage in space (projection and range of corresponding x,y and z coordinates) and time, and which variables it can provide, represented by their CF standard_name (very important!). All this can be seen by printing the reader, or its attributes (e.g. <reader>.variables)

• After intialisation, the reader can be asked (e.g. from an OpenDrift simulation object) to provide a set of variables (e.g. currents, temperatures, zeta etc - again by the CF standard_name's) for a "block" that covers a given set of positions (x,y,z) - which will be the particle positions in the case of an OpenDrift simulation. x and y will here be in the readers native coordinates, but z must be in meters, 0 at the ocean surface (until now), and negative downwards toward seafloor. Thus OpenDrift will take care of transforming x and y to corresponding lon and lat, which is used internally in OpenDrift. However, it is the responsibility of the reader to return data on z-levels in meters, and not in its own native vertical coordinates, as e.g. sigma for ROMS. So the ROMS reader must return blocks/arrays with dimensions ['x', 'y', 'z'] that covers all the requested positions. Therefore roppy.multi_zslice is used to interpolate the sigma-layers to z-layers as hardcoded in the ROMS reader (above).

The above may seems a bit rigid and awkward, but allows the flexibility that any OpenDrift module can be used with any of the Readers. E.g. a drift module can declare that it needs sea_surface_wave_significant_height, and specify in its update method what to do with this variable, without worrying about from where it was obtained. And if a reader declears that it can provide this variable (thanks to consistently using CF standard_name) covering the particles/domain of interest, OpenDrift can request that reader to provide that varible, and perform the interpolation from the returned 3D-block (or 2D, or even 4D in case of ensemble models) onto the particle positions.

Another important detail to mention, is that OpenDrift needs to interpolate in time between "reader blocks" before and after the actual time. And this interpolation requires presently that the z-levels are the same. I.e. it is not possible to interpolate in time between one block with zlevels=[0, -.5, -1, 3] and another block with zlevels=[0, -.5, -1, -2, -3] This is relevant to the question above, as this means that if to allow positive z-values/layers, we must also hardcode some given positive z-values. Note that particles may have any z-value, the above is just the z-layers of the blocks, from which to perform the vertical interpolation. And yes hardcoded z-levels is not very elegant, but I have not found a better solution. Another thing that can be questioned, is whether it is ok to do the vertical interpolation from sigma to z within the ROMS reader, before the later horizontal interpolation in OpenDrift on z-levels.

Sorry for the lengthy explanation, but it might be useful to have a "common operational picture" before proceeding on this task. And I can add some more challenges that now also must be considered:

• Say a particle is at the surface z=1 [masl]. For the next request to the ROMS reader of data for this position, zeta may have lowered to 0.5m. Thus ROMS may be unable to provide z-layers above 0.5m back to OpenDrift. Thus OpenDrift must then extrapolate up to 1m, which is probably straightforward. But say the opposite, a particle is at 0.5 masl, but a the next time step zeta has changed to 1m. Should the particle then be moved up to the new surface, or considered now to be 0.5m below the surface? Again the former is probably the natural, but these things have to be dealt with.
• One important feature of OpenDrift is the robustness of allowing "backup" readers and nesting. Say you are using ROMS native with zeta for a coastal domain, and you have a global ocean model without zeta outside in time and/or space. So coming from the ROMS domain, a particle can have z=-1m and still be at the surface. But entering the outer model domain, the surface will always be at z=0, and the particle will then be considered to be at 1m depth, as this model does not know what the real zeta is (at least not consistent with the inner model).

I believe this latter point is not so easy to deal with. I am very positive to assist in letting OpenDrift support changing elevation and wetting/drying, but several questions like these arise due to the generic nature of OpenDrift. Thus some careful thinking is needed before proceeding.

[knutfrode]

After some more minutes of thinking, and acknowledging that there might be other challenges arising, a pragmatic solution comes to mind. How about:

• z=0 will still always refer to the ocean surface, as a choice of coordinate system / reference. We will then never have particles with z>0.
• whenever zeta is available and not 0, it will be returned as an environment variable (as was is in fact always the case: https://github.com/OpenDrift/opendrift/blob/master/opendrift/readers/reader_ROMS_native.py#L100)
• zeta will then be known for each particle, and its real depth below the sea surface can then be inferred from z, zeta and seafloor depth (H): particle_depth_below_surface = z*(zeta+H)/H
  • I.e. z=0 means depth=0, and z=-H means depth =zeta+H
  • zeta < -H (depth<0) means that we have a dried cell.

Right now it sounds to me that this will omit the new problems, while still allowing to deal properly with the cases of variable zeta combined with shallow ocean.

E.g. to halt a particle on a dried cell will be as simple as: self.elements.moving[self.environment.zeta<=-self.environment.sea_floor_depth_below_sea_level] = 0

moving is an element property indicating whether it should be moved (by currents, turbulence, Stokes...) or not, which can be set back to 1 when water rises: self.elements.moving[self.environment.zeta>-self.environment.sea_floor_depth_below_sea_level] = 1 and is thus different from deactivation.

[kthyng]

@knutfrode I like your pragmatic solution and so far don't see anything wrong with it. I can work try it out from my end today and see how it works.

What functions would use particle_depth_below_surface?

[knutfrode]

Good. In the meantime I am actually convinced that this is by far the best method.

Always defining z=0 as surface makes most of the core things in OpenDrift much easier, without the need to worry about zeta varying in time and space, and being inconsistent between different ocean models:

• when particles are at surface, i.e. when to apply wind drift factor
• applying Stokes drift to the surface, and diminishing with depth downwards
• vertical turbulent eddies that have to be reflected at surface

And then zeta != 0 can always be dealt with as exactly as might be required whenever needed, e.g. for shallow waters and wetting/drying etc. I am not sure whether particle_depth_below_surface might actually be needed, but one should keep in mind that for, say zeta=+1m and H=2m, then z=-1 does actually not mean 1m below surface, but instead z*(zeta+H)/H = -1*(1+3)/2 = -2m below surface. And that seafloor is not where z=-H, but where z=-(zeta+H). Hence z is still depth in m, but scaled.

[kthyng]

Here is how I understand it and see if we agree. Thinking about it in terms of height above seabed because that is straight-forward:

• A drifter depth of z=0 (at the surface) is the same as H + zeta above the seabed
• A drifter depth of z=-1 (1 meter below the time-varying surface) is the same as H + zeta - 1 above the seabed
• A drifter depth at the seabed is the same as H + zeta - H - zeta = 0 above the seabed

I don't understand the scaling you're describing though, could you restate that?

[knutfrode]

Yes, you are right, I made it unnecessary complicated. So

• z=0 is then always ocean surface (possible shifted by zeta)
• z is then the real depth below surface in meters.
• And the only difference that need to be implemented in core OpenDrift is that seafloor should now be where z=-(zeta+H), and not z=-H as presently. But this is only handled in a few centralised places, and should be straightforward. But zeta must now always be available as an environment variable, defaulting to 0.

I believe you mentioned somewhere that you observed some discrepancies between OpenDrift and your own calculations. I guess this could be related to the unfortunate fact that the ROMS reader needs to regrid/interpolate the sigma-layers to regular and hardcoded z-layers, which are then subsequently interpolated vertically within OpenDrift(OceanDrift) e.g. for the vertical turbulent mixing where density profile is needed for each particle. Interpolating each particle/profile directly from the "raw" sigma-profiles would be desireable, but again this is complicated by the fact that OpenDrift is generic and not specific to ROMS.

[kthyng]

Thank you, I should be able to get back to this and other on-going development work on opendrift by mid-next week.

[kthyng]

I am working on another PR that addresses seafloor checks that incorporate zeta and an additional parameter that wetting/drying models have (or at least ROMS does but I think this is commonly necessary) as a minimum wet value — maybe this addresses at least part of your 3rd bullet of the seafloor being z=-(zeta+H). Are there other places that shows up?

As part of this I am adding sea_surface_height to the OceanDrift model as a required variable with a fallback value of 0 and it looks like I need to add it as the same to every model that subclasses OceanDrift — does that sound right? I'll send this in soon and link it into this conversation when I do.

[knutfrode]

Yes, I also believe this means that sea_surface_height (defaulting to 0) must be added to all models subclassing OceanDrift. But this should be fine.

One place where one need to add zeta to H is in the general method interact_with_seafloor https://opendrift.github.io/_modules/opendrift/models/basemodel.html#OpenDriftSimulation.interact_with_seafloor There are probably some other places. However, I think it will not be very critical, as it will not crash, but only lead to an error in the seafloor depth equal to zeta. Thus any other occurances can also be fixed later. So as long as all tests are passing, I think it should be ok for a first merge.

[kthyng]

Yes, interact_with_seafloor is where I've made changes. Now I see I should also modify the vertical_mixing routine in OceanDrift which also checks the seafloor boundary twice and calls interact_with_seafloor. I will do this now.

[kthyng]

@knutfrode I'm writing some tests and had a realization — since z must be <= 0, z_rho also must be <= 0 right? If so I need to rework the formulation I input for sdepth.

[knutfrode]

Sorry for the late feedback. I believe you are right. But this is perhaps now included in #1264 ?

[knutfrode]

Closed by https://github.com/OpenDrift/opendrift/pull/1264