From 6a12d02fcb2fea21d287607007cc5c837899a20f Mon Sep 17 00:00:00 2001 From: johnryantaylor Date: Sun, 18 Oct 2026 22:06:54 +0100 Subject: [PATCH 1/5] Added Dissolving Particles Added DissolvingParticles as a type of active particle. This module can be used to model dissolving OAE feedstock or secondary mineral precipitates. --- docs/make.jl | 4 +- docs/prepare_examples.jl | 3 +- docs/src/appendix/library.md | 6 + .../individuals/dissolving_particles.md | 109 ++++++++ .../src/model_components/individuals/index.md | 2 +- examples/dissolving_particles.jl | 253 ++++++++++++++++++ .../DissolvingParticles.jl | 207 ++++++++++++++ src/Models/Models.jl | 4 + src/OceanBioME.jl | 3 + test/runtests.jl | 1 + test/test_dissolving_particles.jl | 180 +++++++++++++ 11 files changed, 769 insertions(+), 3 deletions(-) create mode 100644 docs/src/model_components/individuals/dissolving_particles.md create mode 100644 examples/dissolving_particles.jl create mode 100644 src/Models/Individuals/DissolvingParticles/DissolvingParticles.jl create mode 100644 test/test_dissolving_particles.jl diff --git a/docs/make.jl b/docs/make.jl index 657b8b15b..93426ad54 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -25,7 +25,8 @@ examples = [ "Baroclinic instability" => "eady", "Model with particles (kelp) interacting with the biogeochemistry" => "kelp", "Data assimilation" => "data_assimilation", - "Ocean alkalinity enhancement" => "oae_experiment" + "Ocean alkalinity enhancement" => "oae_experiment", + "Sinking and dissolving particles" => "dissolving_particles" ] example_pages = [ title => "generated/$(filename).md" for (title, filename) in examples ] @@ -83,6 +84,7 @@ sediments_pages = [ individuals_pages = [ "Overview" => "model_components/individuals/index.md", "Sugar Kelp (Broch and Slagstad 2012 ++)" => "model_components/individuals/slatissima.md", + "Dissolving particles" => "model_components/individuals/dissolving_particles.md", ] component_pages = [ diff --git a/docs/prepare_examples.jl b/docs/prepare_examples.jl index a1dfd3e9d..3d7aded9e 100644 --- a/docs/prepare_examples.jl +++ b/docs/prepare_examples.jl @@ -10,7 +10,8 @@ const EXAMPLES = [ "eady", "kelp", "data_assimilation", - "oae_experiment" + "oae_experiment", + "dissolving_particles" ] for example_name in EXAMPLES diff --git a/docs/src/appendix/library.md b/docs/src/appendix/library.md index 72050ffc0..0804b2531 100644 --- a/docs/src/appendix/library.md +++ b/docs/src/appendix/library.md @@ -96,6 +96,12 @@ Modules = [OceanBioME.Models.PISCESModel.Nitrogen] Modules = [OceanBioME.Models.SugarKelpModel] ``` +### Dissolving particles + +```@autodocs +Modules = [OceanBioME.Models.DissolvingParticlesModel] +``` + ### Carbon Chemistry ```@autodocs diff --git a/docs/src/model_components/individuals/dissolving_particles.md b/docs/src/model_components/individuals/dissolving_particles.md new file mode 100644 index 000000000..537924304 --- /dev/null +++ b/docs/src/model_components/individuals/dissolving_particles.md @@ -0,0 +1,109 @@ +# [Dissolving particles](@id dissolving-particles) + +Dissolving particles are active particles which sink, dissolve, and release material into tracer fields. +They were written for modelling mineral-based ocean alkalinity enhancement (OAE), where particles of a mineral such as +lime (Ca(OH)₂), calcite (CaCO₃), or olivine (Mg₂SiO₄) are added to the ocean and increase the alkalinity as they dissolve, and to model the evolution of secondary mineral precipitates which can form in response to OAE. + +The only property tracked by each particle is its radius, ``r`` (m), which sets how fast the particle dissolves and how fast it sinks. +An example of particles sinking and dissolving in two-dimensional turbulence is given in the [examples](@ref dissolving-particles-example). + +You can set up the particle biogeochemistry with [`DissolvingMineral`](@ref): +```jldoctest +using OceanBioME + +lime = DissolvingMineral() + +# output +DissolvingMineral{Float64} biogeochemistry tracking particle `r`adius +├── density: 2211.0 kg/m³ +├── molar mass: 0.074093 kg/mol +├── released tracers: (Alk = 2.0,) mol/mol +├── dissolution rate: 1.0e-8 m/s +└── sinking speed: StokesSinking{Float64}(2211.0, 1025.0, 0.00108, 9.81) m/s +``` +which can be put into `BiogeochemicalParticles`, or you can directly create the particles with [`DissolvingParticles`](@ref): +```jldoctest +using OceanBioME, Oceananigans + +grid = RectilinearGrid(size = (1, 1, 1), extent = (1, 1, 1)); +particles = DissolvingParticles(10; grid) + +# output +10 BiogeochemicalParticles with DissolvingMineral{Float64} biogeochemistry: +├── fields: (:r,) +└── coupled tracers: (:Alk,) +``` + +The particles can then be added to any biogeochemical model which has the tracers they release, +for example a carbonate system without any biology: +```julia +biogeochemistry = NutrientsPlanktonDetritus(grid; inorganic_carbon = CarbonateSystem(), particles) +``` + +## Model equations + +The radius of each particle shrinks at the dissolution rate, ``k(r)`` (m/s): + +```math +\frac{dr}{dt} = -k(r), +``` + +so the volume of each particle, ``V = 4\pi r^3/3``, decreases at ``4\pi r^2 k(r)``. +When the dissolution is controlled by reactions on the particle surface ``k`` is a constant +(the default), but it can be any function of radius. For example, when dissolution is limited by diffusion +away from the particle surface ``k \propto 1/r``. + +As the particles dissolve they release each tracer at: + +```math +\frac{dX}{dt} = \frac{4\pi r^2 k(r) \nu_X}{V_m}, +``` + +where ``\nu_X`` is the number of moles of tracer ``X`` released per mole of mineral (`released_tracers`), and +``V_m`` is the molar volume of the mineral (`molar_mass / density`). +For example, calcium hydroxide releases two moles of alkalinity per mole (`released_tracers = (; Alk = 2)`), +and calcium carbonate releases two moles of alkalinity and one mole of dissolved inorganic carbon +(`released_tracers = (; Alk = 2, DIC = 1)`). + +The released tracers go into the grid cell nearest to the particle, and each particle can represent many real +particles by setting its `scalefactor`. + +### Sinking + +When the particles are moved by `LagrangianAdvection` (the default) they move with the water, and sink relative to it at `sinking_speed`, +which can be a number or a function of radius. By default they sink at the speed given by Stokes' law +(see [`StokesSinking`](@ref)): + +```math +w = \frac{2(\rho_p - \rho)gr^2}{9\mu}, +``` + +where ``\rho_p`` is the density of the particle, ``\rho`` the density of the water, ``g`` the gravitational acceleration, +and ``\mu`` the dynamic viscosity of the water. This is only valid for small particles (with a Reynolds number, +``2rw\rho/\mu``, less than about one), so a different `sinking_speed` should be given for larger particles. + +Particles which reach the bottom of the domain bounce off it, so they stay within one timestep's sinking of the bottom +and continue to dissolve into the bottom grid cells. + +### Parameter variable names + +| Symbol | Variable name | Units | +|-------------|--------------------|------------------------| +| ``\rho_p`` | `density` | kg / m³ | +| - | `molar_mass` | kg / mol | +| ``\nu_X`` | `released_tracers` | mol / mol | +| ``k`` | `dissolution_rate` | m / s | +| ``w`` | `sinking_speed` | m / s | + +The default `density`, `molar_mass`, and `released_tracers` are for calcium hydroxide. The default dissolution rate is only illustrative +(a 100 μm particle dissolves in about 3 hours) and should be set for the mineral and conditions being modelled. + +## Model conservations + +The total of each released tracer, in the water and in the particles, is conserved. The amount of tracer ``X`` in each particle is: + +```math +X_\text{particle} = \frac{4\pi r^3 \nu_X}{3V_m}. +``` + +Note that because the particle radius and the tracers are stepped forward differently, this is conserved up to a first order timestepping error. When the particles fully dissolve slightly too much tracer is released, about ``0.5k\Delta t / r_0`` of the total with the default `RungeKutta3` timestepper (or ``3k\Delta t / r_0`` with `QuasiAdamsBashforth2`). As a consequence, the timestep should be much shorter than the time it takes the particles to dissolve, ``r_0/k`` to maintain material conservation. diff --git a/docs/src/model_components/individuals/index.md b/docs/src/model_components/individuals/index.md index 2b59427db..f483bc521 100644 --- a/docs/src/model_components/individuals/index.md +++ b/docs/src/model_components/individuals/index.md @@ -1,6 +1,6 @@ # [Individuals](@id individuals) -The effects of individuals can be modelled in OceanBioME. We have implemented this through custom dynamics in the [Lagrangian Particle tracking feature of Oceananigans](https://clima.github.io/OceananigansDocumentation/stable/model_setup/lagrangian_particles/). We have extended these functionalities to make it easier to implement "active" particles which interact with the tracers. We have then implemented a model of [sugar kelp](@ref sugar-kelp) which can be followed as an example of using this functionality. +The effects of individuals can be modelled in OceanBioME. We have implemented this through custom dynamics in the [Lagrangian Particle tracking feature of Oceananigans](https://clima.github.io/OceananigansDocumentation/stable/model_setup/lagrangian_particles/). We have extended these functionalities to make it easier to implement "active" particles which interact with the tracers. We have then implemented a model of [sugar kelp](@ref sugar-kelp) which can be followed as an example of using this functionality, and [dissolving particles](@ref dissolving-particles) which sink and release tracers as they dissolve (e.g. for ocean alkalinity enhancement). To setup particles first create a particle biogeochemistry, e.g.: diff --git a/examples/dissolving_particles.jl b/examples/dissolving_particles.jl new file mode 100644 index 000000000..97b951e9f --- /dev/null +++ b/examples/dissolving_particles.jl @@ -0,0 +1,253 @@ +# # [Sinking and dissolving particles in two-dimensional turbulence](@id dissolving-particles-example) +# +# This example demonstrates how to use [`DissolvingParticles`](@ref). +# This is a simple model of the dissolution of hydrated lime (calcium hydroxide, Ca(OH)₂) are released at the surface of a two-dimensional turbulent flow. +# Each particle sinks and dissolves, releasing alkalinity into the water as it dissolves. +# The dissolution rate and the sinking speed are functions of the particle radius. +# +# This example demonstrates: +# - How to set up [`DissolvingParticles`](@ref) and couple them to a carbonate system +# - How to set up sinking and dissolution rates as functions of particle radius +# - How to check that alkalinity is conserved between the particles and the water + +# ## Install dependencies +# First we check we have the dependencies installed +# ```julia +# using Pkg +# pkg "add OceanBioME, Oceananigans, CairoMakie, JLD2" +# ``` + +# ## Model setup +# First load the required packages +using OceanBioME, Oceananigans, Printf, Random, Statistics +using Oceananigans.Units + +Random.seed!(1234) # Set the random seed for reproducibility +nothing #hide + +# ## Grid +# We use a two-dimensional (x-z) domain which is 64 m wide and periodic in x, and 64 m deep +grid = RectilinearGrid(size = (128, 128), x = (0, 64), z = (-64, 0), + topology = (Periodic, Flat, Bounded)) + +# ## Turbulence +# Start by spinning up some decaying two-dimensional turbulence from random velocities. +# Over the first half an hour the energy moves to larger scales, leaving a field of eddies +# with velocities of a few mm/s, similar to the sinking speeds of the smaller particles +# Create a non-hydrostatic model with WENO advection and no particles or tracers for the spinup. +spinup_model = NonhydrostaticModel(grid; advection = WENO(order = 5)) + +set!(spinup_model, u = (x, z) -> 0.02 * randn(), w = (x, z) -> 0.02 * randn()) + +spinup = Simulation(spinup_model, Δt = 1, stop_time = 30minutes) +conjure_time_step_wizard!(spinup, cfl = 0.5) +run!(spinup) # Run the spinup simulation + +# ## Particles +# Each particle is a sphere of radius ``r``, and [`DissolvingParticles`](@ref) needs two functions of that radius: the `dissolution_rate`, ``k(r)``, which is how fast the radius shrinks (``dr/dt = -k(r)``), and the `sinking_speed`, ``w(r)``, which is how fast the particle sinks relative to the water. +# Either can be a constant number or a function. +# +# Here, we set the dissolution rate so that the mass loss from the particles is proportional to the surface area of the particles. This leads to a constant decrease in the particle radius, ``dr/dt = -k``. Here, set ``k = 10^{-8}`` m/s, so a particle with initial radius ``r_0`` in meters lasts ``r_0 / k`` seconds. +# A different dependence on radius only needs a different function here, for example `dissolution_rate(r) = 5e-13 / r` +# for dissolution limited by diffusion away from the particle surface. +dissolution_rate(r) = 1e-8 # m/s + +# The particles sink at the Stokes' law speed for their size, ``w = 2(\rho_p - \rho)gr^2/9\mu``, so as the particles dissolve they also sink more slowly. +const ρₚ = 2211 # kg/m³, density of calcium hydroxide +const ρ = 1025 # kg/m³, density of seawater +const μ = 1.08e-3 # kg/m/s, dynamic viscosity of seawater +const g = 9.81 # m/s², gravitational acceleration + +sinking_speed(r) = 2 * (ρₚ - ρ) * g * r^2 / (9μ) # m/s, positive downwards + +# The other [`DissolvingMineral`](@ref) parameters (`molar_mass` and `released_tracers`) are left at their defaults, which are for calcium hydroxide. It releases 2 moles of alkalinity per mole dissolved. +# +# We also pass the particle `density`. It does not change how the particles sink, since that is set by `sinking_speed` above, or how fast their radius shrinks, which is set by `dissolution_rate`. +# Instead it sets how much alkalinity is released as the particles dissolve. The particles lose volume at ``4\pi r^2 k(r)``, and the density converts that volume into moles of mineral through the molar volume, ``V_m = M / \rho_p`` (`molar_mass / density`). +# A denser mineral has more moles in each cubic metre, so it releases more alkalinity for the same change in radius. +# If `sinking_speed` were left at its default (by omitting it from mineral_parameters), the density would also set the Stokes' law sinking speed. +# We use the same ``\rho_p`` as in the sinking speed so that the two are consistent. +# The density is also used below to calculate the mass of each particle. +mineral_parameters = (; density = ρₚ, dissolution_rate, sinking_speed) + +# We release 1000 particles with radii between 20 and 80 μm, in a 16 m wide patch at the surface. +n = 1000 + +r₀ = 20e-6 .+ 60e-6 * rand(n) # Set the initial radius of the particles +x₀ = 24 .+ 16 * rand(n) # Set the initial x position of the particles +z₀ = - rand(n) # Set the initial z position of the particles + +particles = DissolvingParticles(n; grid, mineral_parameters) # Create the particles + +# Each computational particle stands for many real particles. +# The number of real particles represented by each computational particle is set by `scalefactor`. +# Here we release 1 kg of lime per metre in the (flat) y direction, and the scale factor is the same for each computational particle. + +mineral = particles.biogeochemistry # Get the mineral parameters + +particle_mass(r) = mineral.density * 4/3 * π * r^3 # Calculate the mass of a particle + +scalefactor = 1 / sum(particle_mass.(r₀)) + +set!(particles, r = r₀, x = x₀, z = z₀, scalefactors = fill(scalefactor, n)) + +particles # Output the particle properties + +# ## Biogeochemistry +# We couple the particles to a carbonate system with dissolved inorganic carbon (`DIC`) and +# alkalinity (`Alk`) without any biology +biogeochemistry = NutrientsPlanktonDetritus(grid; inorganic_carbon = CarbonateSystem(), particles) + +model = NonhydrostaticModel(grid; biogeochemistry, advection = WENO(order = 5)) + +Alk₀ = 2300 # mmol/m³, initial alkalinity + +set!(model, u = spinup_model.velocities.u, w = spinup_model.velocities.w, DIC = 2000, Alk = Alk₀) + +# ## Simulation +# Run the simulation for 2.5 hours, long enough for the largest particles to dissolve, saving the particles and alkalinity field every minute +simulation = Simulation(model, Δt = 1, stop_time = 2.5hours) + +conjure_time_step_wizard!(simulation, cfl = 0.5, max_Δt = 10) + +progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: %s\n", + iteration(sim), prettytime(sim), prettytime(sim.Δt), + prettytime(sim.run_wall_time)) + +add_callback!(simulation, progress_message, TimeInterval(15minutes)) + +filename = "dissolving_particles" + +simulation.output_writers[:fields] = JLD2Writer(model, (; Alk = model.tracers.Alk), + filename = "$filename.jld2", + schedule = TimeInterval(1minute), + overwrite_files = true) + +# Save the particle positions as well as their radii +particle_state(model) = (; x = Array(particles.x), z = Array(particles.z), r = Array(particles.fields.r)) + +simulation.output_writers[:particles] = JLD2Writer(model, (; particles = particle_state), + filename = "$(filename)_particles.jld2", + schedule = TimeInterval(1minute), + overwrite_files = true) + +# ## Run! +run!(simulation) + +# ## Results +# First we load the alkalinity and particle output +Alk = FieldTimeSeries("$filename.jld2", "Alk") + +times = Alk.times +xc, _, zc = nodes(Alk) + +using JLD2 + +file = jldopen("$(filename)_particles.jld2") + +iterations = keys(file["timeseries/t"]) + +x, z, r = ntuple(_ -> zeros(n, length(iterations)), 3) + +for (i, iter) in enumerate(iterations) + x[:, i] = file["timeseries/particles/$iter"].x + z[:, i] = file["timeseries/particles/$iter"].z + r[:, i] = file["timeseries/particles/$iter"].r +end + +close(file) + +# We can now make an animation of the particles sinking and dissolving through the turbulence. +# The particles are shown at their current size (magnified), and disappear when they have dissolved +using CairoMakie + +frame = Observable(1) + +Alkₙ = @lift interior(Alk[$frame], :, 1, :) .- Alk₀ + +remaining = @lift r[:, $frame] .> 0 +xₙ = @lift x[$remaining, $frame] +zₙ = @lift z[$remaining, $frame] +rₙ = @lift r[$remaining, $frame] +sizeₙ = @lift 2e5 * $rₙ + +fig = Figure(size = (700, 600), fontsize = 20) + +title = @lift "t = " * prettytime(times[$frame]) +Label(fig[0, 1:2], title) + +ax = Axis(fig[1, 1]; title = "Alkalinity increase and particles", xlabel = "x (m)", ylabel = "z (m)", + aspect = DataAspect(), limits = ((0, 64), (-64, 0))) +hm = heatmap!(ax, xc, zc, Alkₙ, colormap = Reverse(:deep), colorrange = (0, 80), interpolate = true) +scatter!(ax, xₙ, zₙ; markersize = sizeₙ, color = :white) +Colorbar(fig[1, 2], hm, label = "Alkalinity increase (mmol/m³)") + +record(fig, "$filename.mp4", 1:length(times), framerate = 12) do i + frame[] = i +end +nothing #hide + +# ![](dissolving_particles.mp4) + +# The large particles sink out of the eddies, while the small particles dissolve near the surface +# and are stirred around by the flow. +# +# To see this more clearly, we can compare how deep each particle had sunk when it finished dissolving +# with how deep it would sink in still water. In still water a particle released at the surface +# sinks ``\int_0^{r_0/k} w(r_0 - kt)\,dt = 2(\rho_p - \rho)g r_0^3/27\mu k`` before it dissolves. +dissolved_index = [findfirst(r[p, :] .== 0) for p in 1:n] + +dissolution_depth = [-z[p, dissolved_index[p]] for p in 1:n] + +still_water_depth(r₀) = sinking_speed(r₀) * r₀ / 3dissolution_rate(r₀) + +# We can also look at the horizontally averaged alkalinity increase, and check that the +# alkalinity released into the water is the same as the alkalinity lost from the particles. +# Since the domain is 2D the cell volumes are per metre in the y direction, just like the scalefactors. +using Oceananigans.Operators: volume + +cell_volume = volume(1, 1, 1, grid, Center(), Center(), Center()) + +alkalinity_released = [sum(interior(Alk[i]) .- Alk₀) * cell_volume for i in 1:length(times)] + +molar_volume = mineral.molar_mass / mineral.density # m³/mol +ν = mineral.released_tracers.Alk # mol Alk per mol Ca(OH)₂ + +alkalinity_in_particles = [sum(@. scalefactor * 4/3 * π * r[:, i]^3 / molar_volume * ν * 1000) for i in 1:length(times)] + +total_alkalinity = alkalinity_in_particles[1] + +@info @sprintf("Alkalinity conservation error: %.2e (relative)", + (alkalinity_released[end] + alkalinity_in_particles[end] - total_alkalinity) / total_alkalinity) + +mean_alkalinity_increase = [mean(interior(Alk[i]) .- Alk₀, dims = 1)[1, 1, :] for i in 1:length(times)] + +fig = Figure(size = (1200, 900), fontsize = 20) + +ax1 = Axis(fig[1, 1], xlabel = "Initial radius (μm)", ylabel = "Depth (m)", yreversed = true, + title = "Depth when fully dissolved", limits = (nothing, (0, 64))) +radii = (20:80) * 1e-6 +scatter!(ax1, r₀ * 1e6, dissolution_depth, color = (:black, 0.3), label = "Simulated turbulence") +lines!(ax1, radii * 1e6, still_water_depth.(radii), color = :red, linewidth = 3, label = "Theoretical: still water") +axislegend(ax1, position = :lb) + +ax2 = Axis(fig[1, 3], xlabel = "Time (hours)", ylabel = "Radius (μm)", title = "Particle radius") +for p in 1:10:n + lines!(ax2, times / hour, r[p, :] * 1e6, color = r₀[p], colormap = :viridis, colorrange = (20e-6, 80e-6)) +end + +ax3 = Axis(fig[2, 1], xlabel = "Time (hours)", ylabel = "z (m)", title = "Horizontally averaged alkalinity increase") +hm = heatmap!(ax3, times / hour, zc, stack(mean_alkalinity_increase)', colormap = Reverse(:deep)) +Colorbar(fig[2, 2], hm, label = "mmol/m³") + +ax4 = Axis(fig[2, 3], xlabel = "Time (hours)", ylabel = "Fraction of alkalinity", title = "Alkalinity budget") +lines!(ax4, times / hour, alkalinity_in_particles / total_alkalinity, linewidth = 3, label = "In particles") +lines!(ax4, times / hour, alkalinity_released / total_alkalinity, linewidth = 3, label = "Dissolved") +lines!(ax4, times / hour, (alkalinity_in_particles .+ alkalinity_released) / total_alkalinity, + linewidth = 3, linestyle = :dash, color = :black, label = "Total") +axislegend(ax4, position = :rc) + +save("dissolving_particles.png", fig) +nothing #hide + +# ![](dissolving_particles.png) diff --git a/src/Models/Individuals/DissolvingParticles/DissolvingParticles.jl b/src/Models/Individuals/DissolvingParticles/DissolvingParticles.jl new file mode 100644 index 000000000..27cbf3aac --- /dev/null +++ b/src/Models/Individuals/DissolvingParticles/DissolvingParticles.jl @@ -0,0 +1,207 @@ +""" +This module adds active particles which sink, dissolve, and release material into tracer fields. For example, mineral particles added to the ocean for the purpose of ocean alkalinity enhancement (OAE) or dissolving secondarary mineral precipitate. + +The particles are defined by a radius, r (m), which sets how fast they dissolve and sink. + +Prognostic properties +===================== +* Radius: r (m) + +Tracer dependencies +=================== +None + +Coupled tracers +=============== +Those listed in `released_tracers`, by default alkalinity: Alk (mmol/m³) +""" + +module DissolvingParticlesModel + +export DissolvingMineral, DissolvingParticles, StokesSinking + +using Oceananigans.Grids: AbstractGrid + +using OceanBioME.Particles: BiogeochemicalParticles + +import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers, update_particle_state! +import Oceananigans.Models.LagrangianParticleTracking: particle_w_velocity + +import Base: summary, show + +""" + StokesSinking(FT = Float64; + particle_density, + water_density = 1025, + dynamic_viscosity = 1.08e-3, + gravitational_acceleration = 9.81) + +Function to calculate Stokes' law for the speed (m/s) at which a sphere of radius `r` (m) sinks: + +```math +w = \\frac{2(\\rho_p - \\rho)gr^2}{9\\mu}, +``` + +where ``\\rho_p`` is the `particle_density` (kg/m³), ``\\rho`` the `water_density` (kg/m³), +``\\mu`` the `dynamic_viscosity` (kg/m/s), and ``g`` the `gravitational_acceleration` (m/s²). + +This is only valid when the particle Reynolds number, ``2rw\\rho/\\mu``, is small, +i.e. for particles smaller than about 100 μm. +""" +struct StokesSinking{FT} + particle_density :: FT + water_density :: FT + dynamic_viscosity :: FT + gravitational_acceleration :: FT +end + +StokesSinking(FT = Float64; + particle_density, + water_density = 1025, + dynamic_viscosity = 1.08e-3, + gravitational_acceleration = 9.81) = + StokesSinking{FT}(particle_density, water_density, dynamic_viscosity, gravitational_acceleration) + +@inline function (stokes::StokesSinking)(r) + ρₚ = stokes.particle_density + ρ = stokes.water_density + μ = stokes.dynamic_viscosity + g = stokes.gravitational_acceleration + + return 2 * (ρₚ - ρ) * g * r^2 / (9 * μ) +end + +""" + DissolvingMineral(FT = Float64; + density = 2211, + molar_mass = 74.093e-3, + released_tracers = (; Alk = 2), + dissolution_rate = 1e-8, + sinking_speed = StokesSinking(FT; particle_density = density)) + +Mineral particles which dissolve and release material into tracer fields as they sink. +The only property tracked is the particle radius, `r` (m), which sets how fast the particles dissolve and sink. + +The radius shrinks at the `dissolution_rate`, + +```math +\\frac{dr}{dt} = -k(r), +``` + +so each particle loses volume at a rate of ``4\\pi r^2 k(r)`` and releases tracers at +``4\\pi r^2 k(r) \\nu / V_m``, where ``\\nu`` is the moles of tracer released per mole of mineral (`released_tracers`) and ``V_m`` is the molar volume (`molar_mass / density`). +The released tracers are converted from moles to mmol to match the units of the tracers in OceanBioME. + +Particles also sink relative to the water at `sinking_speed` when they are advected by +`LagrangianAdvection` (the default for `BiogeochemicalParticles`). + +Keyword Arguments +================= + +- `density`: density of the mineral (kg/m³) +- `molar_mass`: molar mass of the mineral (kg/mol) +- `released_tracers`: `NamedTuple` of the moles of each tracer released per mole of mineral dissolved, + e.g. `(; Alk = 2, DIC = 1)` for calcium carbonate +- `dissolution_rate`: the speed (m/s) that the radius shrinks, ``k(r)``, either a number or a function of radius +- `sinking_speed`: the speed (m/s, positive downwards) that the particles sink, either a number or a function of radius + +The default values of `density`, `molar_mass`, and `released_tracers` are for calcium hydroxide (hydrated lime, Ca(OH)₂), and by default the particles sink according to Stokes' law (see [`StokesSinking`](@ref)). +The default `dissolution_rate` of 1e-8 m/s is only illustrative (a 100 μm particle dissolves in about 3 hours), and should be set for the mineral and conditions being modelled. + +Example +======= + +```jldoctest +julia> using OceanBioME + +julia> calcite = DissolvingMineral(density = 2710, molar_mass = 100.09e-3, released_tracers = (; Alk = 2, DIC = 1)) +DissolvingMineral{Float64} biogeochemistry tracking particle `r`adius +├── density: 2710.0 kg/m³ +├── molar mass: 0.10009 kg/mol +├── released tracers: (Alk = 2.0, DIC = 1.0) mol/mol +├── dissolution rate: 1.0e-8 m/s +└── sinking speed: StokesSinking{Float64}(2710.0, 1025.0, 0.00108, 9.81) m/s +``` +""" +struct DissolvingMineral{FT, R, D, W} + density :: FT + molar_mass :: FT + released_tracers :: R + dissolution_rate :: D + sinking_speed :: W +end + +function DissolvingMineral(FT = Float64; + density = 2211, + molar_mass = 74.093e-3, + released_tracers = (; Alk = 2), + dissolution_rate = 1e-8, + sinking_speed = StokesSinking(FT; particle_density = density)) + + released_tracers = map(ν -> convert(FT, ν), released_tracers) + dissolution_rate = maybe_convert(FT, dissolution_rate) + sinking_speed = maybe_convert(FT, sinking_speed) + + return DissolvingMineral(convert(FT, density), convert(FT, molar_mass), released_tracers, dissolution_rate, sinking_speed) +end + +maybe_convert(FT, value::Number) = convert(FT, value) +maybe_convert(FT, func) = func + +# `dissolution_rate` and `sinking_speed` can either be numbers or functions of radius +@inline evaluate(value::Number, r) = value +@inline evaluate(func, r) = func(r) + +""" + DissolvingParticles(n; grid, mineral_parameters = NamedTuple(), kwargs...) + +Sets up `n` dissolving `BiogeochemicalParticles` with default [`DissolvingMineral`](@ref) parameters except those specified in `mineral_parameters`. `kwargs` are passed onto `BiogeochemicalParticles`. +""" +DissolvingParticles(n; grid::AbstractGrid{FT}, mineral_parameters = NamedTuple(), kwargs...) where FT = + BiogeochemicalParticles(n; grid, biogeochemistry = DissolvingMineral(FT; mineral_parameters...), kwargs...) + +@inline required_particle_fields(::DissolvingMineral) = (:r, ) +@inline required_tracers(::DissolvingMineral) = () +@inline coupled_tracers(mineral::DissolvingMineral) = keys(mineral.released_tracers) + +# Volume of mineral dissolving per second (m³/s), which stops once the particle has fully dissolved +@inline dissolution(mineral, r) = ifelse(r > 0, 4π * r^2 * evaluate(mineral.dissolution_rate, r), zero(r)) + +# Change in radius (m/s) +@inline (mineral::DissolvingMineral)(::Val{:r}, t, r) = + - ifelse(r > 0, evaluate(mineral.dissolution_rate, r), zero(r)) + +# Release of each tracer (mmol/s) +@inline function (mineral::DissolvingMineral)(::Val{tracer}, t, r) where tracer + molar_volume = mineral.molar_mass / mineral.density # m³/mol + ν = mineral.released_tracers[tracer] # mol tracer / mol mineral + + return 1000 * ν * dissolution(mineral, r) / molar_volume +end + +@inline sinking_speed(mineral, r) = ifelse(r > 0, evaluate(mineral.sinking_speed, r), zero(r)) + +# Oceananigans' particle advection asks for each particle's vertical velocity, +# so the particles sink relative to the water when they're moved by `LagrangianAdvection` +@inline particle_w_velocity(particles::BiogeochemicalParticles{N, <:DissolvingMineral}, p, w) where N = + w - sinking_speed(particles.biogeochemistry, @inbounds(particles.fields.r[p])) + +# The final timestep of dissolution can overshoot to a negative radius, so we clip it +# (this doesn't change what was released because nothing is released when `r ≤ 0`) +function update_particle_state!(particles::BiogeochemicalParticles{N, <:DissolvingMineral}, model, Δt) where N + r = particles.fields.r + r .= max.(r, 0) + return nothing +end + +summary(::DissolvingMineral{FT}) where FT = "DissolvingMineral{$FT} biogeochemistry" + +show(io::IO, mineral::DissolvingMineral) = + print(io, summary(mineral), " tracking particle `r`adius", "\n", + "├── density: ", mineral.density, " kg/m³", "\n", + "├── molar mass: ", mineral.molar_mass, " kg/mol", "\n", + "├── released tracers: ", mineral.released_tracers, " mol/mol", "\n", + "├── dissolution rate: ", mineral.dissolution_rate, " m/s", "\n", + "└── sinking speed: ", mineral.sinking_speed, " m/s") + +end # module diff --git a/src/Models/Models.jl b/src/Models/Models.jl index d61093ea2..dc227d0db 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -15,6 +15,8 @@ export Oxygen export SugarKelp, SugarKelpParticles, GiantKelp +export DissolvingMineral, DissolvingParticles, StokesSinking + export CarbonChemistry export GasExchange, @@ -31,12 +33,14 @@ include("seawater_density.jl") include("CarbonChemistry/CarbonChemistry.jl") # NPD's ExplicitCalciumCarbonate uses CarbonChemistry, so load it first include("AdvectedPopulations/NutrientsPlanktonDetritus/NutrientsPlanktonDetritus.jl") include("Individuals/SugarKelp/SugarKelp.jl") +include("Individuals/DissolvingParticles/DissolvingParticles.jl") include("GasExchange/GasExchange.jl") include("AdvectedPopulations/PISCES/PISCES.jl") using .SedimentModels using .NutrientsPlanktonDetritusModels using .SugarKelpModel +using .DissolvingParticlesModel using .PISCESModel using .CarbonChemistryModel using .GasExchangeModel diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index a7192f4d3..f0c7a507a 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -26,6 +26,9 @@ export DepthDependantSinkingSpeed, PrescribedLatitude, ModelLatitude, PISCESMode # Macroalgae models export BiogeochemicalParticles, SugarKelp, SugarKelpParticles +# Dissolving particles (e.g. for ocean alkalinity enhancement) +export DissolvingMineral, DissolvingParticles, StokesSinking + # Box model export BoxModel, BoxModelGrid, SpeedyOutput, load_output diff --git a/test/runtests.jl b/test/runtests.jl index eec9ac7ee..259821511 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -5,6 +5,7 @@ include("test_utils.jl") include("test_construction_in_functions.jl") include("test_light.jl") include("test_particles.jl") +include("test_dissolving_particles.jl") #include("test_slatissima.jl") include("test_external_plankton_interface.jl") include("test_NutrientsPlanktonDetritus.jl") diff --git a/test/test_dissolving_particles.jl b/test/test_dissolving_particles.jl new file mode 100644 index 000000000..dc37ca592 --- /dev/null +++ b/test/test_dissolving_particles.jl @@ -0,0 +1,180 @@ +include("dependencies_for_runtests.jl") + +using Oceananigans.Architectures: on_architecture +using OceanBioME.Particles: BiogeochemicalParticles, required_particle_fields, coupled_tracers + +# total tracer in a grid made of 1 m³ cells (mmol) +total_tracer(tracer) = sum(on_architecture(CPU(), interior(tracer))) + +# total of a released tracer held in the particles (mmol) +function particle_content(particles, tracer) + mineral = particles.biogeochemistry + r = on_architecture(CPU(), particles.fields.r) + scalefactors = on_architecture(CPU(), particles.scalefactors) + + molar_volume = mineral.molar_mass / mineral.density + ν = mineral.released_tracers[tracer] + + return sum(@. scalefactors * 4/3 * π * r^3 / molar_volume * ν * 1000) +end + +function dissolving_particle_model(grid, particles) + biogeochemistry = NutrientsPlanktonDetritus(grid; inorganic_carbon = CarbonateSystem(), particles) + + model = NonhydrostaticModel(grid; biogeochemistry, advection = nothing) + + set!(model, DIC = 2000, Alk = 2300) + + return model +end + +@testset "Dissolving particle setup" begin + grid = RectilinearGrid(architecture; size = (2, 2, 2), extent = (2, 2, 2)) + + particles = DissolvingParticles(3; grid) + + @test particles isa BiogeochemicalParticles + @test particles.biogeochemistry isa DissolvingMineral + @test length(particles) == 3 + @test required_particle_fields(particles) == (:r, ) + @test coupled_tracers(particles) == (:Alk, ) + + mineral = DissolvingMineral(Float32; released_tracers = (; Alk = 2, DIC = 1), dissolution_rate = r -> 1e-12 / r) + + @test mineral.density isa Float32 + @test mineral.released_tracers === (Alk = 2f0, DIC = 1f0) + @test mineral.sinking_speed isa StokesSinking{Float32} + @test coupled_tracers(mineral) == (:Alk, :DIC) + + stokes = StokesSinking(; particle_density = 2025, water_density = 1025, + dynamic_viscosity = 1e-3, gravitational_acceleration = 9) + + @test stokes(1e-4) ≈ 2 * 1000 * 9 * 1e-8 / (9 * 1e-3) +end + +@testset "Dissolving particle dissolution and conservation" begin + grid = RectilinearGrid(architecture; size = (2, 2, 2), extent = (2, 2, 2)) + + k = 1e-8 # m/s + r₀ = 1e-4 # m + Δt = 1 + + particles = DissolvingParticles(2; grid, + advection = nothing, + scalefactors = [1e6, 2e6], + mineral_parameters = (; dissolution_rate = k, + released_tracers = (; Alk = 2, DIC = 1))) + + set!(particles, r = r₀, x = [0.5, 1.5], y = 0.5, z = -0.5) + + model = dissolving_particle_model(grid, particles) + + initial_Alk = total_tracer(model.tracers.Alk) + particle_content(particles, :Alk) + initial_DIC = total_tracer(model.tracers.DIC) + particle_content(particles, :DIC) + initial_DIC_in_water = total_tracer(model.tracers.DIC) + initial_Alk_in_water = total_tracer(model.tracers.Alk) + + for _ in 1:100 + time_step!(model, Δt) + end + + # the radius shrinks at the (constant) dissolution rate + @test all(on_architecture(CPU(), particles.fields.r) .≈ r₀ - 100 * k * Δt) + + # released into the cells containing the particles in proportion to their scalefactors + Alk = on_architecture(CPU(), interior(model.tracers.Alk)) + @test Alk[1, 1, 2] > 2300 + @test (Alk[2, 1, 2] - 2300) ≈ 2 * (Alk[1, 1, 2] - 2300) + @test Alk[1, 1, 1] == 2300 + + # released in the right ratio + released_Alk = total_tracer(model.tracers.Alk) - initial_Alk_in_water + released_DIC = total_tracer(model.tracers.DIC) - initial_DIC_in_water + @test released_Alk ≈ 2 * released_DIC + + # conserved up to the timestepping error, which is of order kΔt / r + final_Alk = total_tracer(model.tracers.Alk) + particle_content(particles, :Alk) + final_DIC = total_tracer(model.tracers.DIC) + particle_content(particles, :DIC) + + @test isapprox(final_Alk - initial_Alk, 0; atol = 1e-3 * released_Alk) + @test isapprox(final_DIC - initial_DIC, 0; atol = 1e-3 * released_DIC) +end + +@testset "Dissolving particles dissolve completely" begin + grid = RectilinearGrid(architecture; size = (2, 2, 2), extent = (2, 2, 2)) + + k = 1e-6 # m/s + r₀ = 1e-4 # m + Δt = 0.1 # so the particles dissolve in 1000 timesteps + + particles = DissolvingParticles(1; grid, + advection = nothing, + scalefactors = [1e6], + mineral_parameters = (; dissolution_rate = k)) + + set!(particles, r = r₀, x = 0.5, y = 0.5, z = -0.5) + + model = dissolving_particle_model(grid, particles) + + initial_content = particle_content(particles, :Alk) + + for _ in 1:1100 + time_step!(model, Δt) + end + + # the radius stops at zero rather than overshooting + @test all(on_architecture(CPU(), particles.fields.r) .== 0) + + released_Alk = total_tracer(model.tracers.Alk) - 2300 * 8 + + @test isapprox(released_Alk, initial_content; rtol = 10 * k * Δt / r₀) + + # and nothing more is released once they've dissolved + time_step!(model, Δt) + + @test total_tracer(model.tracers.Alk) - 2300 * 8 == released_Alk +end + +@testset "Dissolving particle sinking" begin + grid = RectilinearGrid(architecture; size = (2, 2, 10), extent = (2, 2, 10)) + + w = 0.01 # m/s + + particles = DissolvingParticles(2; grid, mineral_parameters = (; sinking_speed = w, dissolution_rate = 0)) + + # the second particle has already dissolved so shouldn't sink + set!(particles, r = [1e-4, 0], x = 1, y = 1, z = -1) + + model = dissolving_particle_model(grid, particles) + + for _ in 1:10 + time_step!(model, 10) + end + + z = on_architecture(CPU(), particles.z) + + @test z[1] ≈ -1 - 10 * 10 * w + @test z[2] ≈ -1 + + # particles stay in the domain when they reach the bottom + for _ in 1:100 + time_step!(model, 10) + end + + z = on_architecture(CPU(), particles.z) + + @test -10 <= z[1] < -9.9 + + # by default they sink by Stokes' law + particles = DissolvingParticles(1; grid, mineral_parameters = (; dissolution_rate = 0)) + + set!(particles, r = 5e-5, x = 1, y = 1, z = -1) + + model = dissolving_particle_model(grid, particles) + + time_step!(model, 10) + + z = on_architecture(CPU(), particles.z) + + @test z[1] ≈ -1 - 10 * particles.biogeochemistry.sinking_speed(5e-5) +end From e207a83dc58cf52483b76ade401e2bcce8ee6647 Mon Sep 17 00:00:00 2001 From: johnryantaylor Date: Mon, 19 Oct 2026 22:02:42 +0100 Subject: [PATCH 2/5] Add dissolving_particles to Buildkite examples matrix The docs build expects generated/dissolving_particles.md, but the Buildkite examples step uses its own hardcoded list, so the example was never generated. Co-Authored-By: Claude Opus 5.5 --- .buildkite/pipeline.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index adba1049d..470ef775b 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -100,6 +100,7 @@ steps: - "kelp" - "data_assimilation" - "oae_experiment" + - "dissolving_particles" depends_on: "init-gpu" agents: hostname: "$BUILDKITE_AGENT_META_DATA_HOSTNAME" From f5c35cfad7f3eace0f2113db3ceaf2390a09037c Mon Sep 17 00:00:00 2001 From: johnryantaylor Date: Wed, 28 Oct 2026 20:47:36 +0000 Subject: [PATCH 3/5] review modifications Edits based on Jago's review --- examples/dissolving_particles.jl | 53 +++++------ src/OceanBioME.jl | 1 - test/test_dissolving_particles.jl | 143 ++++++------------------------ 3 files changed, 50 insertions(+), 147 deletions(-) diff --git a/examples/dissolving_particles.jl b/examples/dissolving_particles.jl index 97b951e9f..0bf4fc746 100644 --- a/examples/dissolving_particles.jl +++ b/examples/dissolving_particles.jl @@ -30,19 +30,6 @@ nothing #hide grid = RectilinearGrid(size = (128, 128), x = (0, 64), z = (-64, 0), topology = (Periodic, Flat, Bounded)) -# ## Turbulence -# Start by spinning up some decaying two-dimensional turbulence from random velocities. -# Over the first half an hour the energy moves to larger scales, leaving a field of eddies -# with velocities of a few mm/s, similar to the sinking speeds of the smaller particles -# Create a non-hydrostatic model with WENO advection and no particles or tracers for the spinup. -spinup_model = NonhydrostaticModel(grid; advection = WENO(order = 5)) - -set!(spinup_model, u = (x, z) -> 0.02 * randn(), w = (x, z) -> 0.02 * randn()) - -spinup = Simulation(spinup_model, Δt = 1, stop_time = 30minutes) -conjure_time_step_wizard!(spinup, cfl = 0.5) -run!(spinup) # Run the spinup simulation - # ## Particles # Each particle is a sphere of radius ``r``, and [`DissolvingParticles`](@ref) needs two functions of that radius: the `dissolution_rate`, ``k(r)``, which is how fast the radius shrinks (``dr/dt = -k(r)``), and the `sinking_speed`, ``w(r)``, which is how fast the particle sinks relative to the water. # Either can be a constant number or a function. @@ -85,28 +72,25 @@ particles = DissolvingParticles(n; grid, mineral_parameters) # Create the partic mineral = particles.biogeochemistry # Get the mineral parameters -particle_mass(r) = mineral.density * 4/3 * π * r^3 # Calculate the mass of a particle - -scalefactor = 1 / sum(particle_mass.(r₀)) +scalefactor = 1 / sum(mineral.density * 4/3 * π * r₀.^3) # 1 kg divided by the total particle mass set!(particles, r = r₀, x = x₀, z = z₀, scalefactors = fill(scalefactor, n)) particles # Output the particle properties -# ## Biogeochemistry -# We couple the particles to a carbonate system with dissolved inorganic carbon (`DIC`) and -# alkalinity (`Alk`) without any biology -biogeochemistry = NutrientsPlanktonDetritus(grid; inorganic_carbon = CarbonateSystem(), particles) +# We couple the particles to a carbonate system with dissolved inorganic carbon (`DIC`) and alkalinity (`Alk`) without any biology model = NonhydrostaticModel(grid; biogeochemistry, advection = WENO(order = 5)) +# Set the initial conditions Alk₀ = 2300 # mmol/m³, initial alkalinity +u_ic(x, z) = 0.02 * randn() # m/s +w_ic(x, z) = 0.02 * randn() # m/s -set!(model, u = spinup_model.velocities.u, w = spinup_model.velocities.w, DIC = 2000, Alk = Alk₀) +set!(model, u = u_ic, w = w_ic, DIC = 2000, Alk = Alk₀) -# ## Simulation -# Run the simulation for 2.5 hours, long enough for the largest particles to dissolve, saving the particles and alkalinity field every minute -simulation = Simulation(model, Δt = 1, stop_time = 2.5hours) +# Run the simulation for 3 hours, long enough for the largest particles to dissolve, saving the particles and alkalinity field every minutes +simulation = Simulation(model, Δt = 1, stop_time = 3hours) conjure_time_step_wizard!(simulation, cfl = 0.5, max_Δt = 10) @@ -116,6 +100,15 @@ progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, wall time: add_callback!(simulation, progress_message, TimeInterval(15minutes)) +# After 30 minutes, put the particles, DIC, and alkalinity back to their initial values. +function reset_particles_and_tracers!(sim) + set!(particles, r = r₀, x = x₀, z = z₀, scalefactors = fill(scalefactor, n)) + set!(sim.model, DIC = 2000, Alk = Alk₀, enforce_incompressibility = false) + return nothing +end +add_callback!(simulation, reset_particles_and_tracers!, SpecifiedTimes(30minutes)) + +# Save output filename = "dissolving_particles" simulation.output_writers[:fields] = JLD2Writer(model, (; Alk = model.tracers.Alk), @@ -182,7 +175,10 @@ hm = heatmap!(ax, xc, zc, Alkₙ, colormap = Reverse(:deep), colorrange = (0, 80 scatter!(ax, xₙ, zₙ; markersize = sizeₙ, color = :white) Colorbar(fig[1, 2], hm, label = "Alkalinity increase (mmol/m³)") -record(fig, "$filename.mp4", 1:length(times), framerate = 12) do i +# start the movie after the spinup period has ended +run_frames = findall(times .>= 30minutes) + +record(fig, "$filename.mp4", run_frames, framerate = 12) do i frame[] = i end nothing #hide @@ -203,12 +199,7 @@ still_water_depth(r₀) = sinking_speed(r₀) * r₀ / 3dissolution_rate(r₀) # We can also look at the horizontally averaged alkalinity increase, and check that the # alkalinity released into the water is the same as the alkalinity lost from the particles. -# Since the domain is 2D the cell volumes are per metre in the y direction, just like the scalefactors. -using Oceananigans.Operators: volume - -cell_volume = volume(1, 1, 1, grid, Center(), Center(), Center()) - -alkalinity_released = [sum(interior(Alk[i]) .- Alk₀) * cell_volume for i in 1:length(times)] +alkalinity_released = [Field(Integral(Alk[i] - Alk₀))[1, 1, 1] for i in 1:length(times)] molar_volume = mineral.molar_mass / mineral.density # m³/mol ν = mineral.released_tracers.Alk # mol Alk per mol Ca(OH)₂ diff --git a/src/OceanBioME.jl b/src/OceanBioME.jl index 7320d5154..3545a90e4 100644 --- a/src/OceanBioME.jl +++ b/src/OceanBioME.jl @@ -26,7 +26,6 @@ export DepthDependantSinkingSpeed, PrescribedLatitude, ModelLatitude, PISCESMode # Macroalgae models export BiogeochemicalParticles, SugarKelp, SugarKelpParticles -# Dissolving particles (e.g. for ocean alkalinity enhancement) export DissolvingMineral, DissolvingParticles, StokesSinking # Box model diff --git a/test/test_dissolving_particles.jl b/test/test_dissolving_particles.jl index dc37ca592..9a2bd2399 100644 --- a/test/test_dissolving_particles.jl +++ b/test/test_dissolving_particles.jl @@ -1,22 +1,11 @@ include("dependencies_for_runtests.jl") using Oceananigans.Architectures: on_architecture -using OceanBioME.Particles: BiogeochemicalParticles, required_particle_fields, coupled_tracers # total tracer in a grid made of 1 m³ cells (mmol) total_tracer(tracer) = sum(on_architecture(CPU(), interior(tracer))) -# total of a released tracer held in the particles (mmol) -function particle_content(particles, tracer) - mineral = particles.biogeochemistry - r = on_architecture(CPU(), particles.fields.r) - scalefactors = on_architecture(CPU(), particles.scalefactors) - - molar_volume = mineral.molar_mass / mineral.density - ν = mineral.released_tracers[tracer] - - return sum(@. scalefactors * 4/3 * π * r^3 / molar_volume * ν * 1000) -end +particle_radii(particles) = on_architecture(CPU(), particles.fields.r) function dissolving_particle_model(grid, particles) biogeochemistry = NutrientsPlanktonDetritus(grid; inorganic_carbon = CarbonateSystem(), particles) @@ -28,115 +17,56 @@ function dissolving_particle_model(grid, particles) return model end -@testset "Dissolving particle setup" begin - grid = RectilinearGrid(architecture; size = (2, 2, 2), extent = (2, 2, 2)) - - particles = DissolvingParticles(3; grid) - - @test particles isa BiogeochemicalParticles - @test particles.biogeochemistry isa DissolvingMineral - @test length(particles) == 3 - @test required_particle_fields(particles) == (:r, ) - @test coupled_tracers(particles) == (:Alk, ) - - mineral = DissolvingMineral(Float32; released_tracers = (; Alk = 2, DIC = 1), dissolution_rate = r -> 1e-12 / r) - - @test mineral.density isa Float32 - @test mineral.released_tracers === (Alk = 2f0, DIC = 1f0) - @test mineral.sinking_speed isa StokesSinking{Float32} - @test coupled_tracers(mineral) == (:Alk, :DIC) - - stokes = StokesSinking(; particle_density = 2025, water_density = 1025, - dynamic_viscosity = 1e-3, gravitational_acceleration = 9) - - @test stokes(1e-4) ≈ 2 * 1000 * 9 * 1e-8 / (9 * 1e-3) -end +@testset "Dissolving particle dissolution" begin + grid = RectilinearGrid(architecture; size = (1, 1, 1), extent = (1, 1, 1)) -@testset "Dissolving particle dissolution and conservation" begin - grid = RectilinearGrid(architecture; size = (2, 2, 2), extent = (2, 2, 2)) - - k = 1e-8 # m/s + k = 1e-6 # m/s r₀ = 1e-4 # m - Δt = 1 + Δt = 0.1 # so the particle dissolves in 1000 timesteps + scalefactor = 1e6 - particles = DissolvingParticles(2; grid, + particles = DissolvingParticles(1; grid, advection = nothing, - scalefactors = [1e6, 2e6], + scalefactors = [scalefactor], mineral_parameters = (; dissolution_rate = k, released_tracers = (; Alk = 2, DIC = 1))) - set!(particles, r = r₀, x = [0.5, 1.5], y = 0.5, z = -0.5) + set!(particles, r = r₀, x = 0.5, y = 0.5, z = -0.5) model = dissolving_particle_model(grid, particles) - initial_Alk = total_tracer(model.tracers.Alk) + particle_content(particles, :Alk) - initial_DIC = total_tracer(model.tracers.DIC) + particle_content(particles, :DIC) - initial_DIC_in_water = total_tracer(model.tracers.DIC) - initial_Alk_in_water = total_tracer(model.tracers.Alk) + mineral = particles.biogeochemistry + initial_mineral = scalefactor * 4/3 * π * r₀^3 * mineral.density / mineral.molar_mass * 1000 # mmol - for _ in 1:100 + # the radius shrinks at the dissolution rate + for _ in 1:500 time_step!(model, Δt) end - # the radius shrinks at the (constant) dissolution rate - @test all(on_architecture(CPU(), particles.fields.r) .≈ r₀ - 100 * k * Δt) - - # released into the cells containing the particles in proportion to their scalefactors - Alk = on_architecture(CPU(), interior(model.tracers.Alk)) - @test Alk[1, 1, 2] > 2300 - @test (Alk[2, 1, 2] - 2300) ≈ 2 * (Alk[1, 1, 2] - 2300) - @test Alk[1, 1, 1] == 2300 + @test particle_radii(particles)[1] ≈ r₀ - 500 * k * Δt - # released in the right ratio - released_Alk = total_tracer(model.tracers.Alk) - initial_Alk_in_water - released_DIC = total_tracer(model.tracers.DIC) - initial_DIC_in_water - @test released_Alk ≈ 2 * released_DIC - - # conserved up to the timestepping error, which is of order kΔt / r - final_Alk = total_tracer(model.tracers.Alk) + particle_content(particles, :Alk) - final_DIC = total_tracer(model.tracers.DIC) + particle_content(particles, :DIC) - - @test isapprox(final_Alk - initial_Alk, 0; atol = 1e-3 * released_Alk) - @test isapprox(final_DIC - initial_DIC, 0; atol = 1e-3 * released_DIC) -end - -@testset "Dissolving particles dissolve completely" begin - grid = RectilinearGrid(architecture; size = (2, 2, 2), extent = (2, 2, 2)) - - k = 1e-6 # m/s - r₀ = 1e-4 # m - Δt = 0.1 # so the particles dissolve in 1000 timesteps - - particles = DissolvingParticles(1; grid, - advection = nothing, - scalefactors = [1e6], - mineral_parameters = (; dissolution_rate = k)) - - set!(particles, r = r₀, x = 0.5, y = 0.5, z = -0.5) - - model = dissolving_particle_model(grid, particles) - - initial_content = particle_content(particles, :Alk) - - for _ in 1:1100 + # and stops at zero rather than overshooting + for _ in 1:600 time_step!(model, Δt) end - # the radius stops at zero rather than overshooting - @test all(on_architecture(CPU(), particles.fields.r) .== 0) + @test particle_radii(particles)[1] == 0 - released_Alk = total_tracer(model.tracers.Alk) - 2300 * 8 + # all of the mineral is released, in the right ratio, up to the timestepping error + released_Alk = total_tracer(model.tracers.Alk) - 2300 + released_DIC = total_tracer(model.tracers.DIC) - 2000 - @test isapprox(released_Alk, initial_content; rtol = 10 * k * Δt / r₀) + @test isapprox(released_Alk, 2 * initial_mineral; rtol = 10 * k * Δt / r₀) + @test released_Alk ≈ 2 * released_DIC - # and nothing more is released once they've dissolved + # and nothing more is released once it has dissolved time_step!(model, Δt) - @test total_tracer(model.tracers.Alk) - 2300 * 8 == released_Alk + @test total_tracer(model.tracers.Alk) - 2300 == released_Alk end @testset "Dissolving particle sinking" begin - grid = RectilinearGrid(architecture; size = (2, 2, 10), extent = (2, 2, 10)) + grid = RectilinearGrid(architecture; size = (1, 1, 10), extent = (2, 2, 10)) w = 0.01 # m/s @@ -156,25 +86,8 @@ end @test z[1] ≈ -1 - 10 * 10 * w @test z[2] ≈ -1 - # particles stay in the domain when they reach the bottom - for _ in 1:100 - time_step!(model, 10) - end - - z = on_architecture(CPU(), particles.z) - - @test -10 <= z[1] < -9.9 - - # by default they sink by Stokes' law - particles = DissolvingParticles(1; grid, mineral_parameters = (; dissolution_rate = 0)) - - set!(particles, r = 5e-5, x = 1, y = 1, z = -1) - - model = dissolving_particle_model(grid, particles) - - time_step!(model, 10) - - z = on_architecture(CPU(), particles.z) + stokes = StokesSinking(; particle_density = 2025, water_density = 1025, + dynamic_viscosity = 1e-3, gravitational_acceleration = 9) - @test z[1] ≈ -1 - 10 * particles.biogeochemistry.sinking_speed(5e-5) + @test stokes(1e-4) ≈ 2 * 1000 * 9 * 1e-8 / (9 * 1e-3) end From 52586e0b3353f6a8a8a5849e3c704291208c9a57 Mon Sep 17 00:00:00 2001 From: johnryantaylor Date: Wed, 28 Oct 2026 21:06:42 +0000 Subject: [PATCH 4/5] add field dependence Adds optional required tracers upon which the sinking speed and dissolution rate can depend. E.g. sinking might depend on the water density and dissolution might depend on the saturation state. --- .../individuals/dissolving_particles.md | 18 ++- examples/dissolving_particles.jl | 9 ++ .../DissolvingParticles.jl | 131 +++++++++++++----- test/test_dissolving_particles.jl | 32 ++++- 4 files changed, 150 insertions(+), 40 deletions(-) diff --git a/docs/src/model_components/individuals/dissolving_particles.md b/docs/src/model_components/individuals/dissolving_particles.md index 537924304..23890900e 100644 --- a/docs/src/model_components/individuals/dissolving_particles.md +++ b/docs/src/model_components/individuals/dissolving_particles.md @@ -19,7 +19,8 @@ DissolvingMineral{Float64} biogeochemistry tracking particle `r`adius ├── molar mass: 0.074093 kg/mol ├── released tracers: (Alk = 2.0,) mol/mol ├── dissolution rate: 1.0e-8 m/s -└── sinking speed: StokesSinking{Float64}(2211.0, 1025.0, 0.00108, 9.81) m/s +├── sinking speed: StokesSinking{Float64, Float64}(2211.0, 1025.0, 0.00108, 9.81) m/s +└── required fields: () ``` which can be put into `BiogeochemicalParticles`, or you can directly create the particles with [`DissolvingParticles`](@ref): ```jldoctest @@ -85,6 +86,20 @@ and ``\mu`` the dynamic viscosity of the water. This is only valid for small par Particles which reach the bottom of the domain bounce off it, so they stay within one timestep's sinking of the bottom and continue to dissolve into the bottom grid cells. +### Dependence on model fields + +The dissolution rate and sinking speed can also depend on model fields at the particle position, such as +tracers (e.g. temperature), velocities, or auxiliary fields, by listing their names in `required_fields`. +`dissolution_rate` and `sinking_speed` are then called with the radius followed by the values of those fields. +For example, a dissolution rate which increases with temperature, and Stokes' sinking in water whose density +depends on temperature and salinity: +```julia +mineral = DissolvingMineral(required_fields = (:T, :S), + dissolution_rate = (r, T, S) -> 1e-8 * (1 + 0.05 * (T - 20)), + sinking_speed = StokesSinking(particle_density = 2211, + water_density = (T, S) -> 1025 * (1 - 2e-4 * (T - 10) + 8e-4 * (S - 35)))) +``` + ### Parameter variable names | Symbol | Variable name | Units | @@ -94,6 +109,7 @@ and continue to dissolve into the bottom grid cells. | ``\nu_X`` | `released_tracers` | mol / mol | | ``k`` | `dissolution_rate` | m / s | | ``w`` | `sinking_speed` | m / s | +| - | `required_fields` | - | The default `density`, `molar_mass`, and `released_tracers` are for calcium hydroxide. The default dissolution rate is only illustrative (a 100 μm particle dissolves in about 3 hours) and should be set for the mineral and conditions being modelled. diff --git a/examples/dissolving_particles.jl b/examples/dissolving_particles.jl index 0bf4fc746..a674e0d8b 100644 --- a/examples/dissolving_particles.jl +++ b/examples/dissolving_particles.jl @@ -37,6 +37,10 @@ grid = RectilinearGrid(size = (128, 128), x = (0, 64), z = (-64, 0), # Here, we set the dissolution rate so that the mass loss from the particles is proportional to the surface area of the particles. This leads to a constant decrease in the particle radius, ``dr/dt = -k``. Here, set ``k = 10^{-8}`` m/s, so a particle with initial radius ``r_0`` in meters lasts ``r_0 / k`` seconds. # A different dependence on radius only needs a different function here, for example `dissolution_rate(r) = 5e-13 / r` # for dissolution limited by diffusion away from the particle surface. +# +# The dissolution rate can also depend on model fields at the particle position (tracers, velocities, or auxiliary fields). +# To do this, list the fields in `required_fields` in `mineral_parameters` below, and they are passed after the radius, +# e.g. with `required_fields = (:T, )` we would write `dissolution_rate(r, T) = ...`. dissolution_rate(r) = 1e-8 # m/s # The particles sink at the Stokes' law speed for their size, ``w = 2(\rho_p - \rho)gr^2/9\mu``, so as the particles dissolve they also sink more slowly. @@ -47,6 +51,10 @@ const g = 9.81 # m/s², gravitational acceleration sinking_speed(r) = 2 * (ρₚ - ρ) * g * r^2 / (9μ) # m/s, positive downwards +# `sinking_speed` is also passed the values of any `required_fields`, e.g. `sinking_speed(r, T) = ...`. +# Alternatively, the default [`StokesSinking`](@ref) can take a water density which depends on the +# required fields, e.g. `StokesSinking(particle_density = ρₚ, water_density = (T, S) -> ...)`. + # The other [`DissolvingMineral`](@ref) parameters (`molar_mass` and `released_tracers`) are left at their defaults, which are for calcium hydroxide. It releases 2 moles of alkalinity per mole dissolved. # # We also pass the particle `density`. It does not change how the particles sink, since that is set by `sinking_speed` above, or how fast their radius shrinks, which is set by `dissolution_rate`. @@ -79,6 +87,7 @@ set!(particles, r = r₀, x = x₀, z = z₀, scalefactors = fill(scalefactor, n particles # Output the particle properties # We couple the particles to a carbonate system with dissolved inorganic carbon (`DIC`) and alkalinity (`Alk`) without any biology +biogeochemistry = NutrientsPlanktonDetritus(grid; inorganic_carbon = CarbonateSystem(), particles) model = NonhydrostaticModel(grid; biogeochemistry, advection = WENO(order = 5)) diff --git a/src/Models/Individuals/DissolvingParticles/DissolvingParticles.jl b/src/Models/Individuals/DissolvingParticles/DissolvingParticles.jl index 27cbf3aac..4b999d2c2 100644 --- a/src/Models/Individuals/DissolvingParticles/DissolvingParticles.jl +++ b/src/Models/Individuals/DissolvingParticles/DissolvingParticles.jl @@ -9,7 +9,7 @@ Prognostic properties Tracer dependencies =================== -None +Those listed in `required_fields`, by default none Coupled tracers =============== @@ -20,11 +20,14 @@ module DissolvingParticlesModel export DissolvingMineral, DissolvingParticles, StokesSinking +using Adapt + +using Oceananigans: fields using Oceananigans.Grids: AbstractGrid -using OceanBioME.Particles: BiogeochemicalParticles +using OceanBioME.Particles: BiogeochemicalParticles, LagrangianAdvection, extract_tracer_values -import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers, update_particle_state! +import OceanBioME.Particles: required_particle_fields, required_tracers, coupled_tracers, update_particle_state!, advect_particles! import Oceananigans.Models.LagrangianParticleTracking: particle_w_velocity import Base: summary, show @@ -36,7 +39,8 @@ import Base: summary, show dynamic_viscosity = 1.08e-3, gravitational_acceleration = 9.81) -Function to calculate Stokes' law for the speed (m/s) at which a sphere of radius `r` (m) sinks: +Function to calculate Stokes' law for the speed (m/s) at which a sphere of radius `r` (m) sinks, +called as `stokes(r, field_values...)`: ```math w = \\frac{2(\\rho_p - \\rho)gr^2}{9\\mu}, @@ -45,26 +49,33 @@ w = \\frac{2(\\rho_p - \\rho)gr^2}{9\\mu}, where ``\\rho_p`` is the `particle_density` (kg/m³), ``\\rho`` the `water_density` (kg/m³), ``\\mu`` the `dynamic_viscosity` (kg/m/s), and ``g`` the `gravitational_acceleration` (m/s²). +The `water_density` can be a number or a function of the values of the mineral's `required_fields` +(see [`DissolvingMineral`](@ref)), e.g. `water_density = (T, S) -> ...` with `required_fields = (:T, :S)`. + This is only valid when the particle Reynolds number, ``2rw\\rho/\\mu``, is small, i.e. for particles smaller than about 100 μm. """ -struct StokesSinking{FT} +struct StokesSinking{FT, R} particle_density :: FT - water_density :: FT + water_density :: R dynamic_viscosity :: FT gravitational_acceleration :: FT end -StokesSinking(FT = Float64; - particle_density, - water_density = 1025, - dynamic_viscosity = 1.08e-3, - gravitational_acceleration = 9.81) = - StokesSinking{FT}(particle_density, water_density, dynamic_viscosity, gravitational_acceleration) +function StokesSinking(FT = Float64; + particle_density, + water_density = 1025, + dynamic_viscosity = 1.08e-3, + gravitational_acceleration = 9.81) + + water_density = maybe_convert(FT, water_density) + + return StokesSinking{FT, typeof(water_density)}(particle_density, water_density, dynamic_viscosity, gravitational_acceleration) +end -@inline function (stokes::StokesSinking)(r) +@inline function (stokes::StokesSinking)(r, field_values...) ρₚ = stokes.particle_density - ρ = stokes.water_density + ρ = evaluate(stokes.water_density, field_values...) μ = stokes.dynamic_viscosity g = stokes.gravitational_acceleration @@ -77,7 +88,8 @@ end molar_mass = 74.093e-3, released_tracers = (; Alk = 2), dissolution_rate = 1e-8, - sinking_speed = StokesSinking(FT; particle_density = density)) + sinking_speed = StokesSinking(FT; particle_density = density), + required_fields = ()) Mineral particles which dissolve and release material into tracer fields as they sink. The only property tracked is the particle radius, `r` (m), which sets how fast the particles dissolve and sink. @@ -95,6 +107,10 @@ The released tracers are converted from moles to mmol to match the units of the Particles also sink relative to the water at `sinking_speed` when they are advected by `LagrangianAdvection` (the default for `BiogeochemicalParticles`). +`dissolution_rate` and `sinking_speed` can depend on model fields (tracers, velocities, or auxiliary fields) +at the particle position by listing their names in `required_fields`. The functions are then called with the +radius followed by the field values, e.g. `required_fields = (:T, )` and `dissolution_rate = (r, T) -> ...`. + Keyword Arguments ================= @@ -102,8 +118,9 @@ Keyword Arguments - `molar_mass`: molar mass of the mineral (kg/mol) - `released_tracers`: `NamedTuple` of the moles of each tracer released per mole of mineral dissolved, e.g. `(; Alk = 2, DIC = 1)` for calcium carbonate -- `dissolution_rate`: the speed (m/s) that the radius shrinks, ``k(r)``, either a number or a function of radius -- `sinking_speed`: the speed (m/s, positive downwards) that the particles sink, either a number or a function of radius +- `dissolution_rate`: the speed (m/s) that the radius shrinks, ``k(r)``, either a number or a function of radius (and `required_fields`) +- `sinking_speed`: the speed (m/s, positive downwards) that the particles sink, either a number or a function of radius (and `required_fields`) +- `required_fields`: names of the model fields that `dissolution_rate` and `sinking_speed` depend on The default values of `density`, `molar_mass`, and `released_tracers` are for calcium hydroxide (hydrated lime, Ca(OH)₂), and by default the particles sink according to Stokes' law (see [`StokesSinking`](@ref)). The default `dissolution_rate` of 1e-8 m/s is only illustrative (a 100 μm particle dissolves in about 3 hours), and should be set for the mineral and conditions being modelled. @@ -120,15 +137,21 @@ DissolvingMineral{Float64} biogeochemistry tracking particle `r`adius ├── molar mass: 0.10009 kg/mol ├── released tracers: (Alk = 2.0, DIC = 1.0) mol/mol ├── dissolution rate: 1.0e-8 m/s -└── sinking speed: StokesSinking{Float64}(2710.0, 1025.0, 0.00108, 9.81) m/s +├── sinking speed: StokesSinking{Float64, Float64}(2710.0, 1025.0, 0.00108, 9.81) m/s +└── required fields: () ``` """ -struct DissolvingMineral{FT, R, D, W} +struct DissolvingMineral{FT, R, D, W, F} density :: FT molar_mass :: FT released_tracers :: R dissolution_rate :: D sinking_speed :: W + + # the names of the `required_fields` are a type parameter so that they are known inside kernels + DissolvingMineral(density::FT, molar_mass::FT, released_tracers::R, dissolution_rate::D, sinking_speed::W, + required_fields::Tuple) where {FT, R, D, W} = + new{FT, R, D, W, required_fields}(density, molar_mass, released_tracers, dissolution_rate, sinking_speed) end function DissolvingMineral(FT = Float64; @@ -136,21 +159,23 @@ function DissolvingMineral(FT = Float64; molar_mass = 74.093e-3, released_tracers = (; Alk = 2), dissolution_rate = 1e-8, - sinking_speed = StokesSinking(FT; particle_density = density)) + sinking_speed = StokesSinking(FT; particle_density = density), + required_fields = ()) released_tracers = map(ν -> convert(FT, ν), released_tracers) dissolution_rate = maybe_convert(FT, dissolution_rate) sinking_speed = maybe_convert(FT, sinking_speed) - return DissolvingMineral(convert(FT, density), convert(FT, molar_mass), released_tracers, dissolution_rate, sinking_speed) + return DissolvingMineral(convert(FT, density), convert(FT, molar_mass), released_tracers, dissolution_rate, sinking_speed, + required_fields) end maybe_convert(FT, value::Number) = convert(FT, value) maybe_convert(FT, func) = func -# `dissolution_rate` and `sinking_speed` can either be numbers or functions of radius -@inline evaluate(value::Number, r) = value -@inline evaluate(func, r) = func(r) +# `dissolution_rate`, `sinking_speed`, and `water_density` can either be numbers or functions +@inline evaluate(value::Number, args...) = value +@inline evaluate(func, args...) = func(args...) """ DissolvingParticles(n; grid, mineral_parameters = NamedTuple(), kwargs...) @@ -161,34 +186,65 @@ DissolvingParticles(n; grid::AbstractGrid{FT}, mineral_parameters = NamedTuple() BiogeochemicalParticles(n; grid, biogeochemistry = DissolvingMineral(FT; mineral_parameters...), kwargs...) @inline required_particle_fields(::DissolvingMineral) = (:r, ) -@inline required_tracers(::DissolvingMineral) = () +@inline required_tracers(::DissolvingMineral{<:Any, <:Any, <:Any, <:Any, F}) where F = F @inline coupled_tracers(mineral::DissolvingMineral) = keys(mineral.released_tracers) # Volume of mineral dissolving per second (m³/s), which stops once the particle has fully dissolved -@inline dissolution(mineral, r) = ifelse(r > 0, 4π * r^2 * evaluate(mineral.dissolution_rate, r), zero(r)) +@inline dissolution(mineral, r, field_values...) = + ifelse(r > 0, 4π * r^2 * evaluate(mineral.dissolution_rate, r, field_values...), zero(r)) # Change in radius (m/s) -@inline (mineral::DissolvingMineral)(::Val{:r}, t, r) = - - ifelse(r > 0, evaluate(mineral.dissolution_rate, r), zero(r)) +@inline (mineral::DissolvingMineral)(::Val{:r}, t, r, field_values...) = + - ifelse(r > 0, evaluate(mineral.dissolution_rate, r, field_values...), zero(r)) # Release of each tracer (mmol/s) -@inline function (mineral::DissolvingMineral)(::Val{tracer}, t, r) where tracer +@inline function (mineral::DissolvingMineral)(::Val{tracer}, t, r, field_values...) where tracer molar_volume = mineral.molar_mass / mineral.density # m³/mol ν = mineral.released_tracers[tracer] # mol tracer / mol mineral - return 1000 * ν * dissolution(mineral, r) / molar_volume + return 1000 * ν * dissolution(mineral, r, field_values...) / molar_volume end -@inline sinking_speed(mineral, r) = ifelse(r > 0, evaluate(mineral.sinking_speed, r), zero(r)) +@inline sinking_speed(mineral, r, field_values...) = + ifelse(r > 0, evaluate(mineral.sinking_speed, r, field_values...), zero(r)) -# Oceananigans' particle advection asks for each particle's vertical velocity, -# so the particles sink relative to the water when they're moved by `LagrangianAdvection` -@inline particle_w_velocity(particles::BiogeochemicalParticles{N, <:DissolvingMineral}, p, w) where N = - w - sinking_speed(particles.biogeochemistry, @inbounds(particles.fields.r[p])) +const DissolvingBiogeochemicalParticles = BiogeochemicalParticles{<:Any, <:DissolvingMineral} + +# Oceananigans' particle advection only passes the particles into `particle_w_velocity`, so we +# wrap them with the model fields so that the sinking speed can depend on the `required_fields` +struct ParticlesWithFields{P, G, F, X, Y, Z} + particles :: P + grid :: G + fields :: F + x :: X # read by Oceananigans' particle advection + y :: Y + z :: Z +end + +ParticlesWithFields(particles, grid, fields) = + ParticlesWithFields(particles, grid, fields, particles.x, particles.y, particles.z) + +Adapt.adapt_structure(to, p::ParticlesWithFields) = + ParticlesWithFields(adapt(to, p.particles), adapt(to, p.grid), adapt(to, p.fields)) + +Base.length(p::ParticlesWithFields) = length(p.particles) + +advect_particles!(advection::LagrangianAdvection, particles::DissolvingBiogeochemicalParticles, model, Δt) = + advect_particles!(advection, ParticlesWithFields(particles, model.grid, fields(model)), model, Δt) + +# The particles sink relative to the water when they're moved by `LagrangianAdvection` +@inline function particle_w_velocity(p::ParticlesWithFields, n, w) + particles = p.particles + + r = @inbounds particles.fields.r[n] + field_values = extract_tracer_values(nothing, particles.field_interpolation, particles, p.grid, p.fields, n) + + return w - sinking_speed(particles.biogeochemistry, r, field_values...) +end # The final timestep of dissolution can overshoot to a negative radius, so we clip it # (this doesn't change what was released because nothing is released when `r ≤ 0`) -function update_particle_state!(particles::BiogeochemicalParticles{N, <:DissolvingMineral}, model, Δt) where N +function update_particle_state!(particles::DissolvingBiogeochemicalParticles, model, Δt) r = particles.fields.r r .= max.(r, 0) return nothing @@ -202,6 +258,7 @@ show(io::IO, mineral::DissolvingMineral) = "├── molar mass: ", mineral.molar_mass, " kg/mol", "\n", "├── released tracers: ", mineral.released_tracers, " mol/mol", "\n", "├── dissolution rate: ", mineral.dissolution_rate, " m/s", "\n", - "└── sinking speed: ", mineral.sinking_speed, " m/s") + "├── sinking speed: ", mineral.sinking_speed, " m/s", "\n", + "└── required fields: ", required_tracers(mineral)) end # module diff --git a/test/test_dissolving_particles.jl b/test/test_dissolving_particles.jl index 9a2bd2399..569bbf344 100644 --- a/test/test_dissolving_particles.jl +++ b/test/test_dissolving_particles.jl @@ -7,10 +7,10 @@ total_tracer(tracer) = sum(on_architecture(CPU(), interior(tracer))) particle_radii(particles) = on_architecture(CPU(), particles.fields.r) -function dissolving_particle_model(grid, particles) +function dissolving_particle_model(grid, particles; kwargs...) biogeochemistry = NutrientsPlanktonDetritus(grid; inorganic_carbon = CarbonateSystem(), particles) - model = NonhydrostaticModel(grid; biogeochemistry, advection = nothing) + model = NonhydrostaticModel(grid; biogeochemistry, advection = nothing, kwargs...) set!(model, DIC = 2000, Alk = 2300) @@ -91,3 +91,31 @@ end @test stokes(1e-4) ≈ 2 * 1000 * 9 * 1e-8 / (9 * 1e-3) end + +@testset "Dissolving particles depending on model fields" begin + grid = RectilinearGrid(architecture; size = (1, 1, 10), extent = (2, 2, 10)) + + r₀ = 1e-4 # m + Δt = 1 + + stokes = StokesSinking(; particle_density = 2025, water_density = T -> 1000 + T) + + particles = DissolvingParticles(1; grid, mineral_parameters = (; required_fields = (:T, ), + dissolution_rate = (r, T) -> 1e-9 * T, + sinking_speed = stokes)) + + set!(particles, r = r₀, x = 1, y = 1, z = -1) + + model = dissolving_particle_model(grid, particles; tracers = :T) + + set!(model, T = 10) + + time_step!(model, Δt) + + @test particle_radii(particles)[1] ≈ r₀ - 1e-8 * Δt + + # sinks as if the water density were 1010 kg/m³ (up to the small change in radius during the step) + w = StokesSinking(; particle_density = 2025, water_density = 1010)(r₀) + + @test isapprox(on_architecture(CPU(), particles.z)[1] + 1, - w * Δt; rtol = 1e-3) +end From b45cd5cc00b16fc01b10a89c7d57c215300212fc Mon Sep 17 00:00:00 2001 From: johnryantaylor Date: Sat, 7 Nov 2026 15:36:02 +0000 Subject: [PATCH 5/5] Update runtests.jl updating call to sugar kelp test based on upstream change. --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 259821511..89df59d75 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,7 +6,7 @@ include("test_construction_in_functions.jl") include("test_light.jl") include("test_particles.jl") include("test_dissolving_particles.jl") -#include("test_slatissima.jl") +include("test_sugar_kelp.jl") include("test_external_plankton_interface.jl") include("test_NutrientsPlanktonDetritus.jl") include("test_PISCES.jl")