start to HABs scenario (#64)

* have to run off OpenDrift branch for now
* fixed a bug for running backward in time
* added a new parameter to the ocean model configuration files
* updated tests and docs

Author: kthyng

ci/environment-py3.11.yml

@@ -14,7 +14,7 @@ dependencies:
   - kerchunk ==0.2.7
   - xarray
   - numpy <=1.26.4  # req from opendrift
-  - opendrift ==1.13.0
+  # - opendrift ==1.13.0
   - scipy
   - dask
   - netcdf4
@@ -31,8 +31,8 @@ dependencies:
   - coverage
   - pip
   - pytest
-  # - pip:
+  - pip:
     # - adios_db
-    # - git+https://github.com/OpenDrift/opendrift
+    - git+https://github.com/kthyng/opendrift.git@start_habs
     # - git+https://github.com/fsspec/kerchunk
     # - xroms  # can't be found on conda-forge for CI, but don't need since in runslow tests?

ci/environment-py3.12.yml

@@ -14,7 +14,7 @@ dependencies:
   - kerchunk ==0.2.7
   - xarray
   - numpy <=1.26.4  # req from opendrift
-  - opendrift ==1.13.0
+  # - opendrift ==1.13.0
   - scipy
   - dask
   - netcdf4
@@ -31,7 +31,8 @@ dependencies:
   - coverage
   - pip
   - pytest
-  # - pip:
+  - pip:
+    - git+https://github.com/kthyng/opendrift.git@start_habs
     # - adios_db
     # - git+https://github.com/OpenDrift/opendrift
     # - git+https://github.com/fsspec/kerchunk

ci/environment-py3.13.yml

@@ -14,7 +14,7 @@ dependencies:
   - kerchunk ==0.2.7
   - xarray
   - numpy <=1.26.4  # req from opendrift
-  - opendrift ==1.13.0
+  # - opendrift ==1.13.0
   - scipy
   - dask
   - netcdf4
@@ -31,7 +31,8 @@ dependencies:
   - coverage
   - pip
   - pytest
-  # - pip:
+  - pip:
+    - git+https://github.com/kthyng/opendrift.git@start_habs
     # - adios_db
     # - git+https://github.com/OpenDrift/opendrift
     # - git+https://github.com/fsspec/kerchunk

docs/conf.py

@@ -62,6 +62,11 @@
     "myst_nb",
 ]
 
+myst_enable_extensions = [
+    "dollarmath",  # enables $...$ and $$...$$
+    "amsmath",  # enables \begin{equation}, \text{}, etc.
+]
+
 # for compiling notebooks with mystnb
 # https://docs.readthedocs.io/en/stable/guides/jupyter.html#using-notebooks-in-other-formats
 nb_custom_formats = {

docs/configuration.md

@@ -19,7 +19,7 @@ A handful of `pydantic` BaseModels make up the configuration for PTM. This allow
 
 The main configuration classes are:
 1. `TheManagerConfig`
-1. `OpenDriftConfig` and instances `LarvalFishModelConfig`, `LeewayModelConfig`, `OceanDriftModelConfig`, `OpenOilModelConfig`
+1. `OpenDriftConfig` and instances `LarvalFishModelConfig`, `LeewayModelConfig`, `OceanDriftModelConfig`, `OpenOilModelConfig`, `HarmfulAlgalBloomModelConfig`
 
 Other configuration classes are:
 1. `OceanModelConfig`
@@ -38,35 +38,47 @@ To retrieve the JSON schema for most parameters related to a PTM run, you can ac
 ```{code-cell} ipython3
 import particle_tracking_manager as ptm
 import pprint
+import json
 
-pprint.pprint(ptm.OceanDriftModelConfig.model_json_schema())
+schema = ptm.OceanDriftModelConfig.model_json_schema()
+print(json.dumps(schema, indent=2))
 ```
 
 #### LarvalFishModelConfig
 
 ```{code-cell} ipython3
-pprint.pprint(ptm.LarvalFishModelConfig.model_json_schema())
+schema = ptm.LarvalFishModelConfig.model_json_schema()
+print(json.dumps(schema, indent=2))
 ```
 
 #### LeewayModelConfig
 
 ```{code-cell} ipython3
-pprint.pprint(ptm.LeewayModelConfig.model_json_schema())
+schema = ptm.LeewayModelConfig.model_json_schema()
+print(json.dumps(schema, indent=2))
 ```
 
 #### OpenOilModelConfig
 
 ```{code-cell} ipython3
-pprint.pprint(ptm.OpenOilModelConfig.model_json_schema())
+schema = ptm.OpenOilModelConfig.model_json_schema()
+print(json.dumps(schema, indent=2))
 ```
 
+#### HarmfulAlgalBloomModelConfig
 
-#### TheManageConfig
+```{code-cell} ipython3
+schema = ptm.HarmfulAlgalBloomModelConfig.model_json_schema()
+print(json.dumps(schema, indent=2))
+```
+
+#### TheManagerConfig
 
 You can also examine the schema for `TheManagerConfig` directly, which is a subset of the parameters in the scenario classes (which inherit from `TheManagerConfig`).
 
 ```{code-cell} ipython3
-pprint.pprint(ptm.TheManagerConfig.model_json_schema())
+schema = ptm.TheManagerConfig.model_json_schema()
+print(json.dumps(schema, indent=2))
 ```
 
 
@@ -86,6 +98,7 @@ Though `OpenDrift` has more models available, the currently wrapped `drift_model
 * Leeway: scenario for Search and Rescue of various objects at the surface
 * OpenOil: oil spill scenarios
 * LarvalFish: scenario for fish eggs and larvae that can grow
+* HarmfulAlgalBloom: scenario for modeling harmful algal blooms once they exist to see where they travel or where they came from
 
 Set these with e.g.:
 

docs/index.rst

@@ -25,6 +25,7 @@ To install from PyPI:
    quick_start
    tutorial
    ocean_models
+   scenario_details
    configuration
    plots
 

docs/ocean_models.md

@@ -20,10 +20,11 @@ Here we describe how to user existing ocean model configurations as well as how
 Some ocean models are built into PTM and can be accessed by the input parameter `ocean_model`.
 
 The built-in ocean models are:
-* **NWGOA** (1999–2008) over the Northwest Gulf of Alaska (Danielson, S. L., K. S. Hedstrom, E. Curchitser,	2016. Cook Inlet Model Calculations, Final Report to Bureau of Ocean Energy Management,	M14AC00014,	OCS	Study BOEM 2015-050, University	of Alaska Fairbanks, Fairbanks,	AK,	149 pp.)
+* **CIOFS3** (1999–2024) across Cook Inlet, Alaska, the newest hindcast version of NOAA's CIOFS model. (Thyng, K. M., C. Liu, 2025. Cook Inlet Circulation Modeling - Long-term Hindcast with Improved Freshwater Forcing and Other Attributes, Final Report to the National Oceanic Atmospheric Administration National Centers for Coastal Ocean Science Kasitsna Bay Lab, Axiom Data Science, Anchorage, AK.)
 * **CIOFS** (1999–2022) across Cook Inlet, Alaska, a hindcast version of NOAA's CIOFS model. (Thyng, K. M., C. Liu, M. Feen, E. L. Dobbins, 2023. Cook Inlet Circulation Modeling, Final Report to Oil Spill Recovery Institute, Axiom Data Science, Anchorage, AK.)
 * **CIOFSOP** (mid-2021 through 48 hours from present time) which is the nowcast/forecast version of the CIOFS model. (Shi, L., L. Lanerolle, Y. Chen, D. Cao, R. Patchen, A. Zhang,
 and E. P. Myers, 2020. NOS Cook Inlet Operational Forecast System: Model development and hindcast skill assessment, NOAA Technical Report NOS CS 40, Silver Spring, Maryland, September 2020.)
+* **NWGOA** (1999–2008) over the Northwest Gulf of Alaska (Danielson, S. L., K. S. Hedstrom, E. Curchitser,	2016. Cook Inlet Model Calculations, Final Report to Bureau of Ocean Energy Management,	M14AC00014,	OCS	Study BOEM 2015-050, University	of Alaska Fairbanks, Fairbanks,	AK,	149 pp.)
 
 
 ### Show available ocean models
@@ -41,7 +42,7 @@ ocean_model_registry.all()
 Show each individual ocean_model:
 
 ```{code-cell} ipython3
-ocean_model_registry.show("NWGOA")
+ocean_model_registry.show("CIOFS3")
 ```
 
 ```{code-cell} ipython3
@@ -52,6 +53,10 @@ ocean_model_registry.show("CIOFS")
 ocean_model_registry.show("CIOFSOP")
 ```
 
+```{code-cell} ipython3
+ocean_model_registry.show("NWGOA")
+```
+
 ### Return ocean model object
 
 To instead return the ocean model config object, use `get`:

docs/scenario_details.md

@@ -0,0 +1,86 @@
+# Details about individual drift models
+
+## Harmful Algal Blooms
+
+With the `HarmfulAlgalBloom` model, users can transport an existing bloom forward in time or run backward in time to determine where it originated.
+
+Currently users can model *Pseudo nitzschia* as a near-surface-limited bloom that experiences high mortality when outside viable temperature and salinity ranges. No growth is currently included. This is the start to tracking relative biomass of the particles making up the bloom.
+
+### Biomass and Mortality Framework
+
+The model includes a simple habitat-dependent survival model in which temperature and salinity determine mortality, and biomass evolves according to an exponential decay equation.
+
+#### Environmental Zone Classification
+
+For each particle, the model first evaluates its local temperature and salinity to determine whether environmental conditions fall within a **preferred**, **marginal**, or **lethal** range. This classification is handled by `classify_zone`, which assigns each element to one of:
+
+- `baseline_mortality` (preferred habitat)
+- `medium_mortality` (suboptimal but viable habitat)
+- `high_mortality` (outside survival limits)
+
+Temperature classification uses four thresholds:
+
+- `temperature_pref_min`
+- `temperature_pref_max`
+- `temperature_death_min`
+- `temperature_death_max`
+
+Salinity classification uses an analogous set:
+
+- `salinity_pref_min`
+- `salinity_pref_max`
+- `salinity_death_min`
+- `salinity_death_max`
+
+
+The result is a per-particle assessment of environmental stress.
+
+#### Mortality Rate Selection
+
+The model then assigns each particle a mortality rate according to a tiered decision rule implemented in `choose_mortality_rate`:
+
+1. If **temperature OR salinity** is in a `high_mortality` zone →
+   use `mortality_rate_high`.
+
+2. Else, if **either** variable is in a `medium_mortality` zone
+   (and neither is high) →
+   use `mortality_rate_medium`.
+
+3. Only when **both** temperature and salinity lie in preferred ranges →
+   use `mortality_rate_baseline`.
+
+This "worst-condition wins" approach prevents double-counting stress while ensuring that the most limiting factor controls mortality.
+
+The output is an array `mortality_rates` (units: days$^{-1}$).
+
+#### Biomass Evolution
+
+Biomass is updated using an exponential decay formulation. Growth is not yet implemented, so the net rate is purely negative:
+
+
+This corresponds to integrating:
+
+$$
+\frac{dB}{dt} = (\text{growth_rate} - \text{mortality_rate}) \times B,
+$$
+
+so biomass decays on a timescale determined by the assigned mortality rate. Particles in lethal zones decay fastest; particles in preferred zones decay slowly (not yet implemented). The solution to the equation is
+
+$$
+B(t) = B(0) \exp^{(\text{growth_rate} - \text{mortality_rate})t}
+$$
+
+#### Deactivation of Dead Particles
+
+Once biomass is updated, particles whose biomass falls below the threshold `biomass_dead_threshold` are removed from the simulation (or deactivated). Deactivated particles no longer participate in advection, mixing, or further biological updates, allowing the Lagrangian population to thin naturally in unfavorable environments.
+
+
+
+### Next steps:
+* add growth
+* add medium level and baseline mortality
+* add Alexandrium catenella and Dinophysis spp.
+* add vertical behavior framework - Implement a shared parameterized vertical-movement function (band / diel_band)
+    * Pseudo-nitzschia: shallow band
+    * Alexandrium: diel vertical migration (shallow band by day, deeper band by night)
+    * Dinophysis: mid-depth band (around the pycnocline)

docs/tutorial.md

@@ -53,6 +53,67 @@ m.run_all()
 ```
 
 
+## Harmful Algal Bloom
+
+The goal of this scenario is to examine the transport of an existing harmful algal bloom or determine where an existing bloom originated, respecting temperature and salinity bounds for the species.
+
+### Initialize manager `m`
+
+```{code-cell} ipython3
+m = ptm.OpenDriftModel(drift_model="HarmfulAlgalBloom", lon = -89.8, lat = 29.08,
+                       number=10, steps=40,
+                       ocean_model="TXLA",
+                       start_time="2009-11-19T12:00",
+                       ocean_model_local=False,
+                       species_type="Pseudo_nitzschia",
+                       plots={'spaghetti': {}})
+```
+
+The currently available species are:
+
+```{code-cell} ipython3
+ptm.HarmfulAlgalBloomModelConfig.model_json_schema()["$defs"]["HABSpeciesTypeEnum"]["enum"]
+```
+
+where `custom` is an option to allow the user to input all necessary parameters to represent a species of their choice.
+
+There are parameters available just for the HAB model:
+
+```{code-cell} ipython3
+import json
+print(json.dumps(ptm.models.opendrift.enums.species_types.HABParameters.model_json_schema(), indent=2))
+```
+
+
+The special parameters for each available species are:
+
+```{code-cell} ipython3
+species = ptm.HarmfulAlgalBloomModelConfig.model_json_schema()["$defs"]["HABSpeciesTypeEnum"]["enum"]
+species.remove('custom')
+for specie in species:
+    print(f"{specie}: {ptm.models.opendrift.enums.species_types.SPECIES_HAB_DEFAULTS[specie]}")
+```
+
+The regular parameters that are set for each available species are:
+
+```{code-cell} ipython3
+for specie in species:
+    print(f"{specie}: {ptm.models.opendrift.enums.species_types.SPECIES_HAB_MANAGER_DEFAULTS[specie]}")
+```
+
+The configuration parameters for this simulation are:
+
+```{code-cell} ipython3
+pprint.pprint(m.config.model_dump())
+```
+
+### Run
+
+```{code-cell} ipython3
+m.run_all()
+```
+
+
 ## Leeway (Search and Rescue)
 
 These are simulations of objects that stay at the surface and are transported by both the wind and ocean currents at rates that depend on how much the object sticks up out of and down into the water. The constants to use for those rates have been experimentally determined by the coastguard and are used in this model.

docs/whats_new.md

@@ -1,5 +1,11 @@
 # What's New
 
+## Unreleased
+
+* Added a new model scenario: `HarmfalAlgalBloom`.
+* Fix a bug for running backward in time.
+* Ocean model configuration files now need to include `chunks`
+
 ## v0.13.2 (November 20, 2025)
 
 * If timezone-aware `start_time`, `start_time_end`, or `end_time` is input, the timezone is now removed and it is assumed they are in the same timezone as the model times.