Gridded Environments for Weathering Process

[Gronohibari12]

Hello everyone!
Looking to add weathering processes to a simulation I want to run, and I was wondering if the Water object can be created with Gridded Data. When reading the example script of weathering processes, i see that it's used with static values that don't take into account the spatial distribution of temperature and salinity.
Is there a way to create the Water object with data from netCDF?
I see that the GridTemperature function exists, will that work to create the object, and does something like that exist for salinity?
Thanks a lot!

[jay-hennen]

Unfortunately we only have partial support for weathering with gridded data, specifically with gridded Wind data (environment.environment_objects.GridWind). The Water object doesn't support gridded data, though it has long been on our to-do list.

GridTemperature was a first step in preparing for the implementation of that support.

[ChrisBarker-NOAA]

@jay-hennen -- to be clear GridTemperature exists, but it is not yet possible to connect the weathering code to it?

@Gronohibari12: The good news is that the weathering algorithms are not all that sensitive to temperature an salinity:

Evaporation is pretty sensitive to temperature, but only in a gross way -- e.g. 10C and 25C will result in different evaporation rates, but a couple degrees doesn't matter much.

Salinity only plays a small role in the sedimentation algorithm, and even then, only if there's pretty high sediment load.

So if you pick average values from the time and region of interest, it should work fine (within all the other uncertainties in the process).

That being said -- it would be good to be able to pull that info from the model results -- we'll get there.

[Gronohibari12]

@ChrisBarker-NOAA @jay-hennen Thanks a lot for the replies! So from what I understand even though the Grid Temperature exists its not yet connected to the weathering algorithms. As in my case the temperature values vary from 14C to 23C for the time i want to run the model.
If its not yet connected, i will try to just run it in smaller timeframes.

[ChrisBarker-NOAA]

Yes, that is correct.

14C to 23C is a wide enough range to matter, yes. If you can re-set the water object between runs, then that should work.

Another option -- I haven't tested this, but you should be able to reset the water object mid-run:
Rather than calling Model.full_run(), you can step through the model one time step at a time:

water_obj = gs.Water(15, units={'temperature': 'C'})

# set up model

for step in model: # you can also just call Model.step() until it stops ...
    # step is a dict with holds some information about the model state at least:
    # {'step_num': 1}
    # there may be more in there, depending on the outputters
    new_temp = get_temp_for_timestep(step{'step_num'})
    water_obj.temperature = new_temp

NOTE: we really should output a bit more in step -- like the model's current time.

[Gronohibari12]

@ChrisBarker-NOAA Thanks a lot i will for sure try this out and update on how it went! Just a small question because i probably didnt understand something correctly if i create the loop that updates the temperature and i do a model full run it wont work ?

[ChrisBarker-NOAA]

i do a model full run it wont work ?

You can do a full model run, updating as it goes, but you can't can't call the full_run method, as that will run start to finish without allowing any changes mid run.

my_model.full_run()` 

is essentially just:

for step in my_model:
    pass

but the second lets you examine (and even change) the state of the model as it goes.

There's a bit more (not much) detail here:

https://gnome.orr.noaa.gov/doc/pygnome/scripting/model.html#running-the-model

[Gronohibari12]

Thanks a lot this way it worked just fine!
Is sedimentation implemented yet ?

[Gronohibari12]

Continuing with weatherer issues i get this error :
C:\Users\User.local\share\mamba\envs\gnome\lib\site-packages\gnome\weatherers\spreading.py:534: RuntimeWarning: divide by zero encountered in divide
blob_area_fgv = .5 * (C**2 / (area[s_mask] / vol_frac_le_st[s_mask])) * time_step # make sure area > 0
C:\Users\User.local\share\mamba\envs\gnome\lib\site-packages\gnome\weatherers\spreading.py:536: RuntimeWarning: divide by zero encountered in divide
blob_area_diffusion = ((7. / 6.) * K * ((area[s_mask] / vol_frac_le_st[s_mask]) / K) ** (1. / 7.)) * time_step
C:\Users\User.local\share\mamba\envs\gnome\lib\site-packages\gnome\weatherers\spreading.py:538: RuntimeWarning: invalid value encountered in multiply
new_le_area = area[s_mask] + vol_frac_le_st[s_mask] * (blob_area_fgv + blob_area_diffusion)
Traceback (most recent call last):

File ~.local\share\mamba\envs\gnome\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
exec(code, globals, locals)

File c:\users\user\desktop\python\modeltry.py:115
for step in model:

File ~.local\share\mamba\envs\gnome\lib\site-packages\gnome\model.py:1278 in next
return self.step()

File ~.local\share\mamba\envs\gnome\lib\site-packages\gnome\model.py:1200 in step
self.weather_elements()

File ~.local\share\mamba\envs\gnome\lib\site-packages\gnome\model.py:1065 in weather_elements
w.weather_elements(sc, time_step, model_time)

File ~.local\share\mamba\envs\gnome\lib\site-packages\gnome\weatherers\evaporation.py:207 in weather_elements
self._set_evap_decay_constant(points, model_time, data,

File ~.local\share\mamba\envs\gnome\lib\site-packages\gnome\weatherers\evaporation.py:151 in _set_evap_decay_constant
raise ValueError("Error in Evaporation routine. One of the"

ValueError: Error in Evaporation routine. One of the exponential decay constant is NaN

This is the code that i am currently trying to run is it an issue with the files or maybe the way i average tempature ? (Also will be fixing the imports as adviced in the past)

# Imports
print("Importing GNOME libraries")

import gnome.scripting as gs
from gnome.outputters import NetCDFOutput, ShapeOutput, Renderer
from gnome.movers.random_movers import RandomMover
from gnome.environment import GridWind, GridCurrent, Water, Wind, Waves
from gnome.spills.release import PolygonRelease
from gnome.spills.spill import Spill
from gnome.weatherers import Evaporation, NaturalDispersion, Emulsification
from gnome.spills.substance import Substance, GnomeOil

print("Importing other libraries")
import geopandas as gpd
import datetime
from netCDF4 import Dataset, num2date
import numpy as np

# Temperature function for each timestep
def get_temp_for_timestep(step_num):
    """Retrieve the averaged sea surface temperature for the given model step."""
    if step_num < len(avg_sst_per_hour):
        return avg_sst_per_hour[step_num]
    else:
        return avg_sst_per_hour[-1]  # Return the last available temperature if out of range


print("Setting up data")
# Paths for data files
bna_map = r"C:\Users\User\Desktop\aegeoerasmus\aegeomap.bna"
wind_data = r'C:\Users\User\Desktop\aegeoerasmus\testdata\modwind2.nc'
cur_data = r'C:\Users\User\Desktop\aegeoerasmus\testdata\curdata.nc'
sst_file = r"C:\Users\User\Desktop\aegeoerasmus\testdata\sstdata.nc"
shapefile_path = r"C:\Users\User\Desktop\aegeoerasmus\testdata\aoi1.shp"

print("Initializing the model")
# Model setup
start_time = datetime.datetime(2021, 4, 1, 23, 0, 0)
end_release = datetime.datetime(2021, 4, 9, 23, 0, 0)
time_step = datetime.timedelta(hours=1)
duration = datetime.timedelta(days=6)
model = gs.Model(start_time=start_time, duration=duration, time_step=time_step)

# Load SST data and process average temperatures
print("Loading SST data")
nc_data = Dataset(sst_file, 'r')
sst_data = nc_data.variables['analysed_sst'][:]
time_var = nc_data.variables['time']
times = num2date(time_var[:], units=time_var.units)
sst_data_celsius = sst_data - 273.15  # Convert to Celsius
avg_sst_per_hour = np.mean(sst_data_celsius, axis=(1, 2))  # Average over lat/lon
nc_data.close()

# Initialize water object
water_obj = gs.Water(avg_sst_per_hour[0], units={'temperature': 'C'})

print("Adding environmental data")
# Environmental data
wind = GridWind.from_netCDF(filename=wind_data)
model.environment += wind
model.environment += water_obj

print("Configuring weathering processes")
# Add weatherers
waves = Waves(wind)
model.add_weathering()
model.weatherers += Evaporation(wind=wind, water=water_obj)
model.weatherers += NaturalDispersion(waves=waves, water=water_obj)
model.weatherers += Emulsification(waves=waves)

print("Setting up spill release")
# Load shapefile for polygon spill
gdf = gpd.read_file(shapefile_path)
polygons = gdf.geometry.tolist()
polygon_release = PolygonRelease(
    release_time=start_time,
    end_release_time=end_release,
    polygons=polygons,
    num_per_timestep=80
)



polygon_spill = Spill(
    release=polygon_release,
    amount=40,  # Total amount over entire release period
    units='tonnes',
    substance= GnomeOil(filename=r"C:\Users\User\Desktop\aegeoerasmus\oda_NO00136.json")
)
model.spills += polygon_spill
print("Creating map and movers")
# Map and movers
model.map = gs.MapFromBNA(filename=bna_map, refloat_halflife=-1)
wind_mover = gs.WindMover.from_netCDF(wind_data)
current_mover = gs.CurrentMover.from_netCDF(cur_data)
random_mover = RandomMover(diffusion_coef=100000.0, uncertain_factor=2.0)
model.movers += wind_mover
model.movers += current_mover
model.movers += random_mover

print("Configuring model output")
# Configure output renderer
map_BB = model.map.get_map_bounding_box()
renderer = Renderer(
    map_filename=bna_map,
    output_dir=r"C:\Users\User\Desktop\gnometest",
    output_timestep=time_step,
    map_BB=map_BB
)
model.outputters += renderer
model.outputters += gs.WeatheringOutput(output_dir=r"C:\Users\User\Desktop\gnometest")

print("Starting simulation")
# Run the model with dynamic temperature updates
for step in model:
    step_num = step["step_num"]
    new_temp = get_temp_for_timestep(step_num)
    water_obj.temperature = new_temp  # Update temperature dynamically

print("Simulation complete")

[coconnor8]

Thanks a lot this way it worked just fine! Is sedimentation implemented yet ?

Sedimentation is included in natural dispersion. It is output as a separate category though.

[ChrisBarker-NOAA]

I hope to get some time to. look closer, but a coupel thoughts:

# Add weatherers
waves = Waves(wind)
model.add_weathering()
model.weatherers += Evaporation(wind=wind, water=water_obj)
model.weatherers += NaturalDispersion(waves=waves, water=water_obj)
model.weatherers += Emulsification(waves=waves)

model.add_weathering() adds it all in one step.

And then you are adding them again, so I think you're getting two of each.

you can print out model.weatherers to be sure.

but you should remove the model.add_weathering() call, and sjmiply add them one by one as you have.

Note: the other option is to call only add_weathering, but I"ll need to look more carefully to see what will happen with the water and waves objects ...

Another thing to check is to do the run. without changing the water object, and see if the same error occurs.

There is a chance that the units get confused for. temperature -- internally K is used -- the water object should convert for you, but it's possible that something got lost in the shuffle when you change it after the fact.

[Gronohibari12]

@ChrisBarker-NOAA Thanks a lot for ur reply! I did in fact change the units and removed the extra weatherers but i get the same error. Non of my Data have any NaN values. The problem is probably with how i update the water object as when i use a static value and do a full run everything runs normally.

[ChrisBarker-NOAA]

Sorry -- that was untested -- I'll look into it.

However, I think you've identified another issue -- it seems the weathering isn't working correctly with polygon releases.

We hadn't tested it well, because we usually use it for spills when we are initializing from, e.g. a satellite observation, and the oil is already in an unspecified weathered state, so we don't turn on weathering.

Stay tuned for a fix.

@coconnor8: have I got that right?

[coconnor8]

Yes. Some recent updates to the spreading algorithm moved initializations to the release, but only for point line releases, not polygon releases. It should be an easy fix.

[Gronohibari12]

@ChrisBarker-NOAA Thanks a lot i will be intune fo the fix. U think this is what causes the Error in Evaporation routine. One of the exponential decay constant is NaN? I have found that if i average the SST for every hour the model runs 2 steps and then gives the error if i average 6 hours the model runs 20 steps and then gives the same error.