OceanBioME · jagoosw · May 8, 2026 · May 9, 2026 · May 9, 2026 · May 10, 2026
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "OceanBioME"
 uuid = "a49af516-9db8-4be4-be45-1dad61c5a376"
-version = "0.18.8"
+version = "0.19.0"
 authors = ["Jago Strong-Wright <js2430@damtp.cam.ac.uk> and contributors"]
 
 [deps]
@@ -23,7 +23,7 @@ EnsembleKalmanProcesses = "1, 2"
 GibbsSeaWater = "0.1"
 JLD2 = "0.4, 0.5, 0.6"
 KernelAbstractions = "0.9"
-Oceananigans = "0.105, 0.106, 0.107, 0.108, 0.109, 0.110, 0.111, 0.112"
+Oceananigans = "0.105, 0.106, 0.107, 0.108, 0.109, 0.110, 0.111, 0.112, 0.113"
 OffsetArrays = "1.14"
 SeawaterPolynomials = "0.3, 0.4"
 julia = "1.10"

diff --git a/docs/oceanbiome.bib b/docs/oceanbiome.bib
@@ -357,6 +357,18 @@ @article{Wanninkhof2014
 }
 
 
+@article{Weiss1980,
+	author = {R.F. Weiss and B.A. Price},
+	doi = {https://doi.org/10.1016/0304-4203(80)90024-9},
+	issn = {0304-4203},
+	journal = {Marine Chemistry},
+	number = {5},
+	pages = {347-359},
+	title = {Nitrous oxide solubility in water and seawater},
+	volume = {8},
+	year = {1980}
+}
+
 @article{Weiss1974,
 	author = {R.F. Weiss},
 	doi = {https://doi.org/10.1016/0304-4203(74)90015-2},

diff --git a/docs/src/model_components/air-sea-gas.md b/docs/src/model_components/air-sea-gas.md
@@ -6,7 +6,15 @@ F = k(u_{10}, T)(C_w - C_a),
 ```
 where `k` is the gas transfer velocity.
 
-Our implementation is intended to be generic for any gas, so you can specify `air_concentration`, `water_concentration`, `transfer_velocity`, and `wind_speed` as any function in `GasExchange`, but we also provide constructors and default values for carbon dioxide and oxygen.
+Our implementation is intended to be generic for any gas, so you can specify `air_concentration`, `water_concentration`, and `transfer_velocity` as any function in `GasExchange`, but we also provide constructors and default values for carbon dioxide and oxygen.
+
+The wind speed lives *inside* the transfer velocity rather than alongside it, so that transfer velocities are free to depend on whatever they need to. The default `SchmidtScaledTransferVelocity` scales a `WindSpeedScaledTransferVelocities` — which holds the `wind_speed` and the ``k_{660}(u_{10})`` parameterisation — by the Schmidt number. The carbon dioxide and oxygen constructors still take `wind_speed` directly and build this for you:
+
+```julia
+CO₂_flux = CarbonDioxideGasExchangeBoundaryCondition(; wind_speed = 5)
+```
+
+`wind_speed` (like `air_concentration`) may be a number, a function of `(x, y, t)`, a function of `(i, j, grid, clock, model_fields)` if `discrete_form = true`, a `Field`, or a `FieldTimeSeries`. If you pass a `FieldTimeSeries` which lives on a different grid to the model, also pass `grid = grid` so that it can be wrapped for interpolation onto the model grid.
 
 To setup carbon dioxide and/or oxygen boundary conditions you simply build the condition and then specify it in the model:
 ```@example gasexchange
@@ -40,22 +48,20 @@ where ``c`` is a coefficient (`coeff`) which typically is wind product specific
 
 Currently, the parameters for CO₂ and oxygen are included, but it would be very straightforward to add the parameters given in the original publication for other gases (e.g. inert tracers of other nutrients such as N₂).
 
-### Carbon dioxide partial pressure
+### Carbon dioxide concentration
 
-For most gasses the water concentration `C_w` is simply taken directly from the biogeochemical model or another tracer (in which case `water_concentration` should be set to `TracerConcentration(:tracer_name)`), but for carbon dioxide the fugacity (``fCO_2``) must be derived from the dissolved inorganic carbon (`DIC`) and `Alk`alinity by a `CarbonChemistry` model (please see the docs for [CarbonChemistry](@ref carbon-chemistry)), and used to calculate the partial pressure (``pCO_2``).
+For most gasses the water concentration `C_w` is simply taken directly from the biogeochemical model or another tracer (in which case `water_concentration` should be set to `TracerConcentration(:tracer_name)`), but for carbon dioxide it must be derived from the dissolved inorganic carbon (`DIC`) and `Alk`alinity by a `CarbonChemistry` model (please see the docs for [CarbonChemistry](@ref carbon-chemistry)).
 
-The default parameterisation for the partial pressure (`CarbonDioxideConcentration`) is given by [dickson2007](@citet) and defines the partial pressure to be the mole fraction ``x(CO_2)`` multiplied by the pressure, ``P``, related to the fugacity by:
+The water concentration is the aqueous carbon dioxide concentration in mmol / m³,
 ```math
-fCO_2 = x(CO_2)P\exp\left(\frac{1}{RT}\int_0^P\left(V(CO_2)-\frac{RT}{P'}\right)dP'\right).
+C_w = [CO_2(aq)] = DIC\frac{[H^+]^2}{[H^+]^2 + K_1[H^+] + K_1K_2},
 ```
-The volume (``V``) is related to the gas pressure by the virial expression:
+(`CarbonDioxideConcentration`), and the air concentration is the dry air mole fraction converted onto the same basis by Dalton's law and a solubility,
 ```math
-\frac{PV(CO_2)}{RT}\approx1+\frac{B(x, T)}{V(CO_2)}+\mathcal{O}(V(CO_2)^{-2}),
+C_a = x(CO_2)p_{atm}f_f(T, S)\frac{\rho}{10^3}.
 ```
-and the first virial coefficient ``B`` for carbon dioxide in air can be approximated as:
+The solubility ``f_f`` is the [Weiss1980](@citet) parameterisation (`FF`), which is the solubility ``K_0`` corrected for the water vapour pressure of saturated air and for the non-ideality of the gas phase,
 ```math
-B_{CO_2-\text{air}} \approx B_{CO_2}(T) + 2x(CO_2)\delta_{CO_2-\text{air}}(T),
+f_f = K_0(1 - p_{H_2O})\gamma,
 ```
-where ``\delta`` is the cross virial coefficient.
-
-``B_{CO_2}`` and ``\delta_{CO_2-\text{air}}`` are parameterised by [Weiss1974](@citet) and reccomended in [dickson2007](@citet) as fourth and first order polynomials respectively.
+and is therefore the right quantity to multiply a *dry air* mole fraction by. The transfer velocity is a bare piston velocity (its `solubility` is `UnitSolubility`), and the atmospheric pressure enters the flux exactly once, on the air side.
diff --git a/docs/src/model_components/carbon-chemistry.md b/docs/src/model_components/carbon-chemistry.md
@@ -328,8 +328,11 @@ The chemical system described above has a large number of equilibrium constants,
 By default, this model parameterises them based on the "best practice" guidelines of [dickson2007](@citet).
 These parameterisations are:
 ```@example carbon-chem
-using OceanBioME.Models.CarbonChemistryModel: K0, K1, K2, KB, KW, KS, KF, KP1, KP2, KP3, KSi, KSP_calcite, KSP_aragonite
-K0() # Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)
+using OceanBioME.Models.CarbonChemistryModel: FF, K0, K1, K2, KB, KW, KS, KF, KP1, KP2, KP3, KSi, KSP_calcite, KSP_aragonite
+K0() # Weiss (1974, Mar. Chem., 2, 203-215)
+```
+```@example carbon-chem
+FF() # Weiss & Price (1980, Mar. Chem., 8, 347-359; Eq 13 with table 6 values)
 ```
 ```@example carbon-chem
 K1() # Millero (1995, Geochim. Cosmochim. Acta, 59, 664)

diff --git a/docs/src/model_components/light.md b/docs/src/model_components/light.md
@@ -81,4 +81,24 @@ using OceanBioME
 light_attenuation = PrescribedAttenuationPAR(grid, surface_PAR; attenuation = 0.1)
 ```
 
-where `surface_PAR` may be a constant or a function `f(x, y, t)`, and `attenuation` may be a constant or a function `f(x, y, z, t)`. `PrescribedAttenuationPAR` also accepts the `interface_field` keyword described in [Recording PAR at cell faces](@ref).
+where `surface_PAR` may be a constant or a function `f(x, y, t)`, and `attenuation` may be a constant or a function `f(x, y, z, t)`. `PrescribedAttenuationPAR` also accepts the `interface_field` keyword described in [Recording PAR at cell faces](@ref).
+
+## Diagnosing the surface PAR from shortwave radiation
+
+All of the models above take a `surface_PAR`, which is usually a constant or a function of horizontal position and time. When you have a shortwave radiation flux instead — for example from an atmospheric forcing dataset, or from a coupled model — `PARFromShortwave` wraps it and takes a fixed fraction of it:
+
+```math
+PAR_0 = f_{PAR}Q_{sw},
+```
+
+where ``f_{PAR}`` is `photosynthetic_fraction_of_shortwave`, ``0.43`` by default. It is used in place of any other `surface_PAR`:
+
+```julia
+using OceanBioME
+
+light_attenuation = TwoBandPhotosyntheticallyActiveRadiation(grid, PARFromShortwave(shortwave))
+```
+
+where `shortwave` may be a constant, a function, a `Field`, or a `FieldTimeSeries`.
+
+When OceanBioME is coupled to [NumericalEarth](https://github.com/NumericalEarth/NumericalEarth.jl), `PARFromShortwave(grid)` builds the surface field for you, and the coupled model writes the ocean's penetrating shortwave radiation (i.e. after reflection and any sea ice blocking) into it every time step.
diff --git a/src/Light/2band.jl b/src/Light/2band.jl
@@ -98,15 +98,7 @@ function TwoBandPhotosyntheticallyActiveRadiation(grid::AbstractGrid{FT}, surfac
     chlorophyll_blue_exponent = convert(FT, chlorophyll_blue_exponent)
     pigment_ratio = convert(FT, pigment_ratio)
 
-    boundary_condition_kwargs = surface_PAR isa Function ? (; parameters, discrete_form) : NamedTuple()
-
-    field = CenterField(grid; boundary_conditions =
-                            regularize_field_boundary_conditions(
-                                FieldBoundaryConditions(top = ValueBoundaryCondition(surface_PAR; boundary_condition_kwargs...)), grid, :PAR))
-
-    # wrap surface_PAR to make it work with the `getbc` interface
-    surface_PAR = materialize_condition(surface_PAR, parameters, discrete_form, ()) 
-    surface_PAR = regularize_boundary_condition(surface_PAR, grid, (Center(), Center(), Center()), 3, RightBoundary, nothing)
+    field, surface_PAR = PAR_field(grid, surface_PAR, parameters, discrete_form)
 
     return TwoBandPhotosyntheticallyActiveRadiation(water_red_attenuation,
                                                     water_blue_attenuation,

diff --git a/src/Light/Light.jl b/src/Light/Light.jl
@@ -6,14 +6,15 @@ module Light
 export TwoBandPhotosyntheticallyActiveRadiation, 
        PrescribedPhotosyntheticallyActiveRadiation, 
        MultiBandPhotosyntheticallyActiveRadiation,
-       PrescribedAttenuationPAR
+       PrescribedAttenuationPAR,
+       PARFromShortwave
 
 using Adapt
 
 using KernelAbstractions, Oceananigans.Units
 using Oceananigans.Architectures: device, architecture, on_architecture
 using Oceananigans.Utils: launch!
-using Oceananigans: Center, Face, fields
+using Oceananigans: Oceananigans, Center, Face, fields
 using Oceananigans.Grids: node, znodes, znode, AbstractGrid
 using Oceananigans.Fields: CenterField, TracerFields, location
 using Oceananigans.BoundaryConditions: fill_halo_regions!, 
@@ -37,6 +38,22 @@ import Base: show, summary
 import Oceananigans.Biogeochemistry: biogeochemical_auxiliary_fields, update_biogeochemical_state!, required_biogeochemical_auxiliary_fields
 import Oceananigans.BoundaryConditions: _fill_top_halo!
 
+function PAR_field(grid, surface_PAR, parameters, discrete_form)
+    boundary_condition_kwargs = surface_PAR isa Function ? (; parameters, discrete_form) : NamedTuple()
+
+    boundary_conditions =
+       regularize_field_boundary_conditions(
+           FieldBoundaryConditions(top = ValueBoundaryCondition(surface_PAR; boundary_condition_kwargs...)), grid, :PAR)
+
+    field = CenterField(grid; boundary_conditions)
+
+    # wrap surface_PAR to make it work with the `getbc` interface
+    surface_PAR = materialize_condition(surface_PAR, parameters, discrete_form, ())
+    surface_PAR = regularize_boundary_condition(surface_PAR, grid, (Center(), Center(), Center()), 3, RightBoundary, nothing)
+
+    return field, surface_PAR
+end
+
 include("abstract_light.jl")
 include("2band.jl")
 include("multi_band.jl")
@@ -45,4 +62,6 @@ include("prescribed_attenuation.jl")
 
 include("compute_euphotic_depth.jl")
 
+include("PAR_from_shortwave.jl")
+
 end
diff --git a/src/Light/PAR_from_shortwave.jl b/src/Light/PAR_from_shortwave.jl
@@ -0,0 +1,47 @@
+import Oceananigans.BoundaryConditions: getbc
+
+"""
+    PARFromShortwave(surface_shortwave; photosynthetic_fraction_of_shortwave = 0.43)
+
+A surface photosynthetically active radiation (`PAR`) which is diagnosed as a fixed fraction of the
+surface downwelling shortwave radiation, i.e.
+
+```math
+PAR_0 = f_{PAR} Q_{sw},
+```
+
+where ``f_{PAR}`` is `photosynthetic_fraction_of_shortwave`.
+
+`surface_shortwave` may be anything the boundary condition `getbc` interface accepts (a number, a
+`Field`, a `FieldTimeSeries`, ...), and is passed as the `surface_PAR` of any light attenuation model,
+e.g. `TwoBandPhotosyntheticallyActiveRadiation(grid, PARFromShortwave(Qsw))`.
+
+When coupled to NumericalEarth, `PARFromShortwave(grid)` builds the surface field for you and the
+coupled model writes the ocean's penetrating shortwave into it each step.
+"""
+struct PARFromShortwave{SS, FT}
+                       surface_shortwave :: SS
+    photosynthetic_fraction_of_shortwave :: FT
+end
+
+Adapt.adapt_structure(to, ad::PARFromShortwave) =
+    PARFromShortwave(adapt(to, ad.surface_shortwave),
+                     adapt(to, ad.photosynthetic_fraction_of_shortwave))
+
+PARFromShortwave(surface_shortwave;
+                 photosynthetic_fraction_of_shortwave = 0.43) =
+    PARFromShortwave(surface_shortwave, photosynthetic_fraction_of_shortwave)
+
+@inline Oceananigans.BoundaryConditions.getbc(light::PARFromShortwave, i, j, args...) =
+    @inbounds light.photosynthetic_fraction_of_shortwave * getbc(light.surface_shortwave, i, j, args...)
+
+summary(::PARFromShortwave) = "PARFromShortwave"
+
+function show(io::IO, light::PARFromShortwave)
+    msg = "PARFromShortwave\n"
+    msg *= "└── photosynthetic_fraction_of_shortwave: $(light.photosynthetic_fraction_of_shortwave)\n"
+
+    print(io, msg)
+
+    return nothing
+end
diff --git a/src/Light/abstract_light.jl b/src/Light/abstract_light.jl
@@ -1,6 +1,10 @@
 using Oceananigans.Operators: Δzᵃᵃᶜ
 
-abstract type AbstractSingleBandExponentialLightAttenuation{BA, IN, CA, SP} end # bands, interface, cell average, and surface light
+abstract type AbstractPhotosyntheticallyActiveRadiation{SP} end
+
+surface_PAR(par::AbstractPhotosyntheticallyActiveRadiation) = par.surface_PAR
+
+abstract type AbstractSingleBandExponentialLightAttenuation{BA, IN, CA, SP} <: AbstractPhotosyntheticallyActiveRadiation{SP} end # bands, interface, cell average, and surface light
 
 const AbstractLight{BA, IN, CA, SP} = AbstractSingleBandExponentialLightAttenuation{BA, IN, CA, SP}
 

diff --git a/src/Light/multi_band.jl b/src/Light/multi_band.jl
@@ -23,7 +23,7 @@ and χ(i) the chlorophyll attenuation coefficient.
 When the fields are called with `biogeochemical_auxiliary_fields` an additional field named `PAR`
 is also returned which is a sum of the bands.
 """
-struct MultiBandPhotosyntheticallyActiveRadiation{T, F, FN, K, E, C, SPAR, SPARD}
+struct MultiBandPhotosyntheticallyActiveRadiation{T, F, FN, K, E, C, SPAR, SPARD} <: AbstractPhotosyntheticallyActiveRadiation{SPAR}
                                 total :: T
                                fields :: F
                           field_names :: FN

diff --git a/src/Light/prescribed.jl b/src/Light/prescribed.jl
@@ -22,7 +22,7 @@ fields which are user specified, e.g. they may be `FunctionField`s or
 fields which will be returned in `biogeochemical_auxiliary_fields`, if only
 one field is present the field will be named `PAR`.
 """
-struct PrescribedPhotosyntheticallyActiveRadiation{F}
+struct PrescribedPhotosyntheticallyActiveRadiation{F} <: AbstractPhotosyntheticallyActiveRadiation{nothing}
     fields :: F
 
     function PrescribedPhotosyntheticallyActiveRadiation(fields)

diff --git a/src/Light/prescribed_attenuation.jl b/src/Light/prescribed_attenuation.jl
@@ -16,7 +16,7 @@ Alias for [`PrescribedAttenuationPhotosyntheticallyActiveRadiation`](@ref).
 const PrescribedAttenuationPAR = PrescribedAttenuationPhotosyntheticallyActiveRadiation
 
 @inline attenuation(i, j, k, grid, la::PrescribedAttenuationPAR, clock, chlorophyll) =
-    la.attenuation(i, j, k, grid, clock, nothing)
+    la.attenuation(i, j, k, grid, clock, chlorophyll)
 
 """
     PrescribedAttenuationPAR(grid, surface_PAR;
@@ -62,16 +62,7 @@ function PrescribedAttenuationPhotosyntheticallyActiveRadiation(grid, surface_PA
                                                                 attenuation_discrete_form = false,
                                                                 interface_field = nothing)
 
-    boundary_condition_kwargs = surface_PAR isa Function ? (; parameters = surface_parameters, discrete_form = surface_discrete_form) : NamedTuple()
-
-    field = CenterField(grid; 
-                        boundary_conditions =
-                            regularize_field_boundary_conditions(
-                                FieldBoundaryConditions(top = ValueBoundaryCondition(surface_PAR; boundary_condition_kwargs...)), grid, :PAR
-                            ))
-
-    surface_PAR = materialize_condition(surface_PAR, surface_parameters, surface_discrete_form, ()) 
-    surface_PAR = regularize_boundary_condition(surface_PAR, grid, (Center(), Center(), Center()), 3, RightBoundary, nothing)
+    field, surface_PAR = PAR_field(grid, surface_PAR, surface_parameters, surface_discrete_form)
 
     if attenuation isa Number
         attenuation = Forcing(ConstantField(attenuation))

diff --git a/...dvectedPopulations/NutrientsPlanktonDetritus/InorganicCarbon/abstract_inorganic_carbon.jl b/...dvectedPopulations/NutrientsPlanktonDetritus/InorganicCarbon/abstract_inorganic_carbon.jl
@@ -1,4 +1,4 @@
-abstract type AbstractInorganicCarbon end
+abstract type AbstractInorganicCarbon{N} end
 
 const NPD_AIC{FT} = NutrientsPlanktonDetritus{FT, <:Any, <:Any, <:Any, <:AbstractInorganicCarbon}
 

diff --git a/...vectedPopulations/NutrientsPlanktonDetritus/InorganicCarbon/explicit_calcium_carbonate.jl b/...vectedPopulations/NutrientsPlanktonDetritus/InorganicCarbon/explicit_calcium_carbonate.jl
@@ -142,7 +142,7 @@ Keyword Arguments
 - `open_bottom`: whether `CaCO₃` can sink out of the bottom of the domain
 - `carbon_chemistry`: the [`CarbonChemistry`](@ref) used to compute the calcium carbonate saturation `Ω`
 """
-struct ExplicitCalciumCarbonate{N, R, CC, SV, SS} <: AbstractInorganicCarbon
+struct ExplicitCalciumCarbonate{N, R, CC, SV, SS} <: AbstractInorganicCarbon{N}
           calcium_carbonate_dissolution_rate :: R
       calcium_carbonate_dissolution_exponent :: R
         calcium_carbonate_precipitation_rate :: R

diff --git a/...vectedPopulations/NutrientsPlanktonDetritus/InorganicCarbon/implicit_calcium_carbonate.jl b/...vectedPopulations/NutrientsPlanktonDetritus/InorganicCarbon/implicit_calcium_carbonate.jl
@@ -10,7 +10,7 @@ Passing `replicates > 1` manifests `replicates` independent copies of the carbon
 (`DIC1`, `Alk1`, `DIC2`, …), which is useful for ensemble or perturbation experiments; each replicate
 evolves with the same tendency as the base `DIC`/`Alk`.
 """
-struct CarbonateSystem{N} <: AbstractInorganicCarbon end
+struct CarbonateSystem{N} <: AbstractInorganicCarbon{N} end
 
 function CarbonateSystem(replicates = 1)
     manifest_carbonate_replicates!(replicates)

diff --git a/src/Models/CarbonChemistry/equilibrium_constants.jl b/src/Models/CarbonChemistry/equilibrium_constants.jl
@@ -83,6 +83,49 @@ summary(::IO, ::K0) = string("Solubility constant")
 show(io::IO, k0::K0) = print(io, "Solubility constant\n",
     "    ln(k₀/k°) = $(k0.constant) + $(k0.inverse_T) / T + $(k0.log_T) (log(T) - log(100)) + $(k0.T²) T² + ($(k0.S) + $(k0.ST) T + $(k0.ST²) T²)S")
 
+"""
+    FF(; constant = -162.8301,
+         inverse_T =  218.2968 * 100,
+         log_T =  90.9241,
+         T² = -1.47696 / 100^2,
+         S =  0.025695,
+         ST = -0.025225 / 100,
+         ST² =  0.0049867 / 100^2)
+
+Parameterisation for the carbon dioxide solubility used to convert a *dry air* mole
+fraction into a dissolved concentration,
+
+    ff = K₀ (1 - pH₂O) γ,
+
+i.e. the solubility `K0` corrected for the water vapour pressure of saturated air and for
+the non-ideality of the gas phase. It is therefore the quantity to multiply a dry air mole
+fraction by (as opposed to `K0`, which multiplies a fugacity), and it has the same units as
+`K0` (mol / kg / atm).
+
+Default values from Weiss, R.F. and Price, B.A. (1980, Mar. Chem., 8, 347-359), equation 13
+with the table 6 values.
+"""
+@kwdef struct FF{FT}
+     constant :: FT = -162.8301
+    inverse_T :: FT =  218.2968 * 100
+        log_T :: FT =  90.9241
+           T² :: FT = -1.47696 / 100^2
+            S :: FT =  0.025695
+           ST :: FT = -0.025225 / 100
+          ST² :: FT =  0.0049867 / 100^2
+end
+
+@inline (c::FF)(T::FT, S; P = nothing) where FT =
+    exp(c.constant
+        + c.inverse_T / T
+        + c.log_T * (log(T) - log(convert(FT, 100)))
+        + c.T² * T^convert(FT, 2)
+        + (c.S + c.ST * T + c.ST² * T^convert(FT, 2)) * S)
+
+summary(::IO, ::FF) = string("Dry air solubility constant")
+show(io::IO, ff::FF) = print(io, "Dry air solubility constant\n",
+    "    ln(ff/k°) = $(ff.constant) + $(ff.inverse_T) / T + $(ff.log_T) (log(T) - log(100)) + $(ff.T²) T² + ($(ff.S) + $(ff.ST) T + $(ff.ST²) T²)S")
+
 """
     K1(FT = Float64;
        constant =  61.2172,