Help configuring Parcels for 3D simulations with MARS3D data

[YunaCreach]

Question

Dear all,

I am working on a project studying larval dispersion to understand connectivity in coastal environments. I am using hydrodynamic data from the MARS3D model (https://mars3d.ifremer.fr/docs/doc.basic.summary.html), a coastal hydrodynamic model that uses sigma vertical coordinates.

The depth z is not directly provided in the MARS3D outputs. My goal is to simulate 3D particle trajectories, releasing particles at a specific depth. However, I am unsure if I need to explicitly calculate z or if there is a more efficient way to handle this in Parcels. Additionally, the hc parameter varies spatially, which complicates using existing Parcels examples (e.g., CROCO).

Currently, the script below does not produce correct simulations. I am reaching out for guidance on how to properly configure Parcels for 3D particle tracking using the outputs of MARS3D

Don't hesitate to ask if anything needs further explanation. Thanks!

# Open the NetCDF file
file = 'C:/CDD/2013/resu_parcels_2013_1h_V11_2013_06.nc'
ds=xr.open_dataset(file)
print(ds.variables)

# Define variables and file paths
filenames = {
    "U": file,
    "V": file,
    "W": file,
    "H": file,
    "Zeta": file,
    "Cs_w": file,
}
variables = {
    "U": "UZ",
    "V": "VZ",
    "W": "w",
    "H": "H0",
    "Zeta": "XE",
    "Cs_w": "SIG_W",
}
dimensions = {
    "U": {"lon": "longitude", "lat": "latitude", "depth": "level_w", "time": "time"},
    "V": {"lon": "longitude", "lat": "latitude", "depth": "level_w", "time": "time"},
    "W": {"lon": "longitude", "lat": "latitude", "depth": "level_w", "time": "time"},
    "H": {"lon": "longitude", "lat": "latitude"},
    "Zeta": {"lon": "longitude", "lat": "latitude", "time": "time"},
    "Cs_w": {"depth": "level_w"},
}

# Create the FieldSet
fieldset = parcels.FieldSet.from_croco(
    file,
    variables,
    dimensions,
    hc=xr.open_dataset(
        file
    ).hc.values,
    allow_time_extrapolation=True,
    mesh="flat",
)

# Define a custom particle deletion function
def DeleteParticle(particle, fieldset, time):
    if particle.state == parcels.StatusCode.ErrorOutOfBounds:
        particle.delete()

# Load spawning locations
spawning_locs = pd.read_csv('U:/CDD/ORM/Parcels/Sites_ponte.csv', sep=';', encoding='latin1')
lat = spawning_locs.lat.tolist()
lon = spawning_locs.lon.tolist()

lon_particles, lat_particles, release_dates = [], [], []
radius = 0.05  # Radius in degrees
dates = [datetime(2013, 6, 3)]

# Generate particle release points
for center_lon, center_lat in zip(lon, lat):
    for date in dates:
        for _ in range(100):  # Particles per site
            angle = np.random.uniform(0, 2 * np.pi)
            r = np.sqrt(np.random.uniform(0, 1)) * radius
            lon_particles.append(center_lon + r * np.cos(angle))
            lat_particles.append(center_lat + r * np.sin(angle))
            release_dates.append(date)

# Create ParticleSet
pset = parcels.ParticleSet(
    fieldset=fieldset,
    pclass=parcels.JITParticle,
    lon=lon_particles,
    lat=lat_particles,
    depth=-20,
    time=release_dates
)

# Run the particle simulation
output_file = pset.ParticleFile(name="Particles.zarr", outputdt=timedelta(hours=1))

pset.execute(
    [parcels.AdvectionRK4_3D, DeleteParticle],
    runtime=timedelta(days=6),
    dt=timedelta(minutes=5),
    output_file=output_file,
)

[erikvansebille]

Thanks for reaching out @YunaCreach. Unfortunately, I've never used the MARS3D model myself, so have no experience with its data. Perhaps others can help, but I can imagine that providing a bit more information can steer the discussion. Perhaps also answer the following questions?

1. In what form is hc given? Is it a full 2D field?
2. What do you mean with "the script below does not produce correct simulations"? What is not correct?

[YunaCreach]

Hi Erik,

Thank you for your reply! To answer your questions:

Form of hc: It is a 2D field. As soon as you move away from the coastline, the value is exactly 9. However, near the coastline, the values can be negative (e.g., -4). I tried fixing the value at 9, but it did not work either.

Issue with the script: The script doesn’t work correctly because the depth values in the fieldset range between 0 and 1, which I suppose is linked to the level-w (sigma).

Let me know if further clarification would help. I’d greatly appreciate any suggestions you or others might have!