🖥️ CROCO tutorial#
This tutorial will show how to run a 3D simulation with output from the CROCO model.
Example setup#
We start with loading the relevant modules and the data.
import matplotlib.pyplot as plt
import numpy as np
import polars as pl
import parcels
import parcels.tutorial
ds_fields = parcels.tutorial.open_dataset("CROCOidealized_data/data")
The simulation will be a very simple, idealised flow: a purely zonal flow over a sloping bottom. This flow (which is somewhat unrealistic of course) nicely showcases that particles stay on their initial depth levels, even though the sigma-layers slope down.
This flow has been created by first running the example from the Shelf front example on the CROCO website. Then, we took the restart file are manually replaced all the u-velocities with 1 m/s and all the v-velocities with 0 m/s. This way we get a purely zonal flow. We then started a new simulation from the restart file, and CROCO then automatically calculated the w velocities to match the new zonal field. We saved the time=0 snapshot from this new run and use it below.
Now we create a FieldSet object using the parcels.convert.croco_to_sgrid() function to first create an S-Grid compliant datatset, and then use that in parcels.FieldSet.from_sgrid_conventions() to create the FieldSet.
Note
In the code below, we use the w velocity field for vertical velocity. However, it is unclear whether this is always the right choice. CROCO (and ROMS) also output an omega field, which may be more appropriate to use. The idealised simulation below only works when using w, though. In other simulations, it is recommended to test whether omega provides more realistic results. See https://github.com/Parcels-code/Parcels/discussions/1728 for more information.
fields = {
"U": ds_fields["u"],
"V": ds_fields["v"],
"W": ds_fields["w"],
"h": ds_fields["h"],
"zeta": ds_fields["zeta"],
"Cs_w": ds_fields["Cs_w"],
}
ds_fset = parcels.convert.croco_to_sgrid(fields=fields, coords=ds_fields)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
# Add the critical depth (`hc`) as context to the fieldset
fieldset.add_context("hc", ds_fields.hc.item())
# Convert the FieldSet to windowed arrays for better performance
fieldset = fieldset.to_windowed_arrays()
INFO: Converted time coordinate from float to timedelta based on units.
DEBUG: Dataset contains 'time' dimension but no 'axis' attribute. Setting 'axis' attribute to 'T'.
Now we can use this Fieldset to advect particles as we would normally do. Note that the particle depths should be provided in (negative) meters, not in sigma-coordinates.
X, Z = np.meshgrid(
[40e3, 80e3, 120e3],
[100, -10, -130, -250, -400, -850, -1400, -1550],
)
Y = np.ones(X.size) * 98000
def DeleteParticle(particles, fieldset):
any_error = particles.state >= 50
particles[any_error].state = parcels.StatusCode.Delete
pset = parcels.ParticleSet(fieldset=fieldset, pclass=parcels.Particle, x=X, y=Y, z=Z)
outputfile = parcels.ParticleFile(
path="croco_particles3D.parquet",
outputdt=np.timedelta64(5000, "s"),
)
pset.execute(
[parcels.kernels.AdvectionRK2_3D_CROCO, DeleteParticle],
runtime=np.timedelta64(50_000, "s"),
dt=np.timedelta64(100, "s"),
output_file=outputfile,
)
Now we plot the particle trajectories below. Note that the particles stay on their initial depth levels, even though the sigma-layers slope down. Also note that particles released above the surface (where depth >0) or below the bathymetry are not advected (due to the DeleteParticle kernel).
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
df = pl.read_parquet("croco_particles3D.parquet")
ax.plot(X / 1e3, Z, "k.", label="Initial positions")
for traj in df.partition_by("particle_id", maintain_order=True):
ax.plot(traj["x"] / 1e3, traj["z"], ".-")
for z in ds_fields.s_w.values:
ax.plot(
fieldset.h.grid.lon[0, :] / 1e3,
fieldset.h.data[0, 0, 25, :] * z,
"k",
linewidth=0.5,
)
ax.set_xlabel("X [km]")
ax.set_xlim(30, 170)
ax.set_ylabel("Depth [m]")
ax.set_title("Particles in idealized CROCO velocity field using AdvectionRK2_3D_CROCO")
plt.tight_layout()
plt.show()
A CROCO simulation with no vertical velocities#
It may be insightful to compare this 3D run with the AdvectionRK2_3D_CROCO kernel with a run where the vertical velocity (W) is set to zero. In that case, the particles will not stay on their initial depth levels but instead follow sigma-layers.
fieldset_noW = parcels.FieldSet.from_sgrid_conventions(ds_fset)
fieldset_noW.W.data[:] = 0.0
fieldset_noW.add_context("hc", ds_fields.hc.item())
fieldset_noW = fieldset_noW.to_windowed_arrays()
X, Z = np.meshgrid(
[40e3, 80e3, 120e3],
[100, -10, -130, -250, -400, -850, -1400, -1550],
)
Y = np.ones(X.size) * 100
pset_noW = parcels.ParticleSet(
fieldset=fieldset_noW, pclass=parcels.Particle, x=X, y=Y, z=Z
)
outputfile = parcels.ParticleFile(
path="croco_particles_noW.parquet", outputdt=np.timedelta64(5000, "s")
)
pset_noW.execute(
[parcels.kernels.AdvectionRK2_3D_CROCO, DeleteParticle],
runtime=np.timedelta64(50_000, "s"),
dt=np.timedelta64(100, "s"),
output_file=outputfile,
)
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
df = pl.read_parquet("croco_particles_noW.parquet")
ax.plot(X / 1e3, Z, "k.", label="Initial positions")
for traj in df.partition_by("particle_id", maintain_order=True):
ax.plot(traj["x"] / 1e3, traj["z"], ".-")
for z in ds_fields.s_w.values:
ax.plot(
fieldset.h.grid.lon[0, :] / 1e3,
fieldset.h.data[0, 0, 25, :] * z,
"k",
linewidth=0.5,
)
ax.set_xlabel("X [km]")
ax.set_xlim(30, 170)
ax.set_ylabel("Depth [m]")
ax.set_title("Particles in idealized CROCO velocity field with W=0")
plt.tight_layout()
plt.show()
The algorithms used#
When using Croco data, Parcels needs to convert the particles.z to sigma-coordinates, before doing any interpolation. This is done with the convert_z_to_sigma_croco() function. Interpolating onto a Field is then done like:
def SampleTempCroco(particles, fieldset):
"""Sample temperature field on a CROCO sigma grid by first converting z to sigma levels."""
sigma = convert_z_to_sigma_croco(fieldset, particles.t, particles.z, particles.y, particles.x, particles)
particles.temp = fieldset.T[particles.t, sigma, particles.y, particles.x, particles]
For Advection, you will need to use the AdvectionRK2_3D_CROCO kernel, which works slightly different from the normal 3D advection kernel because it converts the vertical velocity in sigma-units. The conversion from z to sigma is done at every time step, using the convert_z_to_sigma_croco() function.
In particular, the following algorithm is used (note that the RK2 version is slightly more complex than this Euler-Forward version, but the idea is identical). Also note that the vertical velocity is linearly interpolated here, which gives much better results than the default C-grid interpolation.
sigma = particles.z / fieldset.h[particles.t, 0, particles.y, particles.x]
sigma_levels = convert_z_to_sigma_croco(fieldset, particles.t, particles.z, particles.y, particles.x, particles)
(u, v) = fieldset.UV[particles.t, sigma_levels, particles.y, particles.x, particles]
w = fieldset.W[particles.t, sigma_levels, particles.y, particles.x, particles]
# scaling the w with the sigma level of the particle
w_sigma = w * sigma / fieldset.h[particles.t, 0, particles.y, particles.x, particles]
x_new = particles.x + u*particles.dt
y_new = particles.y + v*particles.dt
# calculating new sigma level
sigma_new = sigma + w_sigma*particles.dt
# converting back from sigma to z, at _new_ location
z_new = sigma_new * fieldset.h[particles.t, 0, y_new, x_new, particles]