diff --git a/core/src/include/ModelComponent.hpp b/core/src/include/ModelComponent.hpp
index 73e944c94..578b492f4 100644
--- a/core/src/include/ModelComponent.hpp
+++ b/core/src/include/ModelComponent.hpp
@@ -1,7 +1,7 @@
 /*!
  * @file ModelComponent.hpp
  *
- * @date 1 Jul 2024
+ * @date 30 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  * @author Einar Ólason <einar.olason@nersc.no>
  */
@@ -71,6 +71,8 @@ namespace Protected {
         = "SLAB_QDW"; // Slab ocean temperature nudging heat flux, W m⁻²
     inline constexpr TextTag SLAB_FDW
         = "SLAB_FDW"; // Slab ocean salinity nudging water flux, kg s⁻¹ m⁻²
+    inline constexpr TextTag FWFLUX = "FWFLUX"; // Fresh-water flux at the ocean surface, kg m⁻²
+    inline constexpr TextTag SFLUX = "SFLUX"; // Salt flux at the ocean surface, kg m⁻²
 }
 
 namespace Shared {
@@ -88,6 +90,9 @@ namespace Shared {
     inline constexpr TextTag DQIA_DT
         = "DQIA_DT"; // Derivative of Qᵢₐ w.r.t. ice surface temperature  W m⁻² K⁻¹
     inline constexpr TextTag Q_PEN_SW = "Q_PEN_SW"; // Penetrating shortwave flux W m⁻²
+    inline constexpr TextTag Q_SW_OW = "Q_SW_OW"; // Net shortwave flux at ocean surface W m⁻²
+    inline constexpr TextTag Q_SW_BASE
+        = "Q_SW_BASE"; // Net shortwave flux at the base of the ice W m⁻²
     // Mass fluxes
     inline constexpr TextTag HSNOW_MELT = "HSNOW_MELT"; // Thickness of snow that melted, m
     // Atmospheric conditions
diff --git a/physics/src/SlabOcean.cpp b/physics/src/SlabOcean.cpp
index bed9d2936..815019c99 100644
--- a/physics/src/SlabOcean.cpp
+++ b/physics/src/SlabOcean.cpp
@@ -1,7 +1,7 @@
 /*!
  * @file SlabOcean.cpp
  *
- * @date 7 Sep 2023
+ * @date 30 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -75,33 +75,29 @@ std::unordered_set<std::string> SlabOcean::hFields() const { return { sstSlabNam
 void SlabOcean::update(const TimestepTime& tst)
 {
     double dt = tst.step.seconds();
+
     // Slab SST update
     qdw = (sstExt - sst) * cpml / relaxationTimeT;
     HField qioMean = qio * cice; // cice at start of TS, not updated
     HField qowMean = qow * (1 - cice); // 1- cice = open water fraction
     sstSlab = sst - dt * (qioMean + qowMean - qdw) / cpml;
+
     // Slab SSS update
     HField arealDensity = cpml / Water::cp; // density times depth, or cpml divided by cp
     // This is simplified compared to the finiteelement.cpp calculation
     // Fdw = delS * mld * physical::rhow /(timeS*M_sss[i] - ddt*delS) where delS = sssSlab - sssExt
     fdw = (1 - sssExt / sss) * arealDensity / relaxationTimeS;
-    // ice volume change, both laterally and vertically
-    HField deltaIceVol = newIce + deltaHice * cice;
-    // change in snow volume due to melting (should be < 0)
-    HField meltSnowVol = deltaSmelt * cice;
-    // Mass per unit area after all the changes in water volume
-    HField denominator
-        = arealDensity - deltaIceVol * Ice::rho - meltSnowVol * Ice::rhoSnow - (emp - fdw) * dt;
-    // Clamp the denominator to be at least 1 m deep, i.e. at least Water::rho kg m⁻²
-    denominator.clampAbove(Water::rho);
-    // Effective ice salinity is always less than or equal to the SSS
-    HField effectiveIceSal = sss;
-    effectiveIceSal.clampBelow(Ice::s);
-    sssSlab = sss
-        + ((sss - effectiveIceSal) * Ice::rho * deltaIceVol // Change due to ice changes
-              + sss * meltSnowVol
-              + (emp - fdw) * dt) // snow melt, precipitation and nudging fluxes.
-            / denominator;
+
+    /* We use the "virtual salinity flux" approach:
+     *   \partial S / \partial t = F_{fw} S / H,
+     * where S is the ocean salinity, F_{fw} is the freshwater flux and H the mixed-layer depth. In
+     * the presence of sea ice F_{fw} S must be replaced by (S-S_i)I + S(P-E+M_s), with I and S
+     * the fresh-water flux due to ice and snow melt, respectively, and P and E, as precipitation
+     * and evaporation, respectively. Using the salt (F_s) and fresh-water fluxes (F_{fw})
+     * calculated by the IOceanBoundary class, we can write the change in ocean salinity as
+     *   \partial S / \partial t = ( S F_{fw} - 10^3 F_s ) / H
+     */
+    sssSlab = sss + (sss * (fwFlux + fdw) - 1e3 * sFlux) * dt / arealDensity;
 }
 
 } /* namespace Nextsim */
diff --git a/physics/src/include/SlabOcean.hpp b/physics/src/include/SlabOcean.hpp
index 804ea059c..20074efde 100644
--- a/physics/src/include/SlabOcean.hpp
+++ b/physics/src/include/SlabOcean.hpp
@@ -1,7 +1,7 @@
 /*!
  * @file SlabOcean.hpp
  *
- * @date 7 Sep 2023
+ * @date 30 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -38,9 +38,8 @@ class SlabOcean : public ModelComponent, public Configured<SlabOcean> {
         , cice(getStore())
         , qio(getStore())
         , qow(getStore())
-        , newIce(getStore())
-        , deltaHice(getStore())
-        , deltaSmelt(getStore())
+        , fwFlux(getStore())
+        , sFlux(getStore())
     {
     }
 
@@ -82,9 +81,8 @@ class SlabOcean : public ModelComponent, public Configured<SlabOcean> {
     ModelArrayRef<Shared::Q_IO, RW> qio;
     ModelArrayRef<Shared::Q_OW, RW> qow;
     // TODO ModelArrayRef to assimilation flux
-    ModelArrayRef<Shared::NEW_ICE, RW> newIce;
-    ModelArrayRef<Shared::DELTA_HICE, RW> deltaHice;
-    ModelArrayRef<Shared::HSNOW_MELT, RW> deltaSmelt;
+    ModelArrayRef<Protected::FWFLUX, RO> fwFlux;
+    ModelArrayRef<Protected::SFLUX, RO> sFlux;
 
     static const std::string sstSlabName;
     static const std::string sssSlabName;
diff --git a/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp b/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp
index bd859d26e..80df0441d 100644
--- a/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp
+++ b/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp
@@ -1,7 +1,7 @@
 /*!
  * @file FiniteElementFluxes.cpp
  *
- * @date Apr 29, 2022
+ * @date 02 Sep 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -56,7 +56,6 @@ void FiniteElementFluxes::setData(const ModelState::DataMap& ms)
     evap.resize();
     Q_lh_ow.resize();
     Q_sh_ow.resize();
-    Q_sw_ow.resize();
     Q_lw_ow.resize();
     Q_lh_ia.resize();
     Q_sh_ia.resize();
@@ -137,22 +136,29 @@ void FiniteElementFluxes::calculateIce(size_t i, const TimestepTime& tst)
     Q_lh_ia[i] = subl[i] * latentHeatIce(tice.zIndexAndLayer(i, 0));
     double dmdot_dT = dragIce_t * rho_air[i] * v_air[i] * dshice_dT[i];
     double dQlh_dT = latentHeatIce(tice.zIndexAndLayer(i, 0)) * dmdot_dT;
+
     // Sensible heat flux
     Q_sh_ia[i]
         = dragIce_t * rho_air[i] * cp_air[i] * v_air[i] * (tice.zIndexAndLayer(i, 0) - t_air[i]);
     double dQsh_dT = dragIce_t * rho_air[i] * cp_air[i] * v_air[i];
+
     // Shortwave flux
     double albedoValue, i0;
     std::tie(albedoValue, i0)
         = iIceAlbedoImpl->surfaceShortWaveBalance(tice.zIndexAndLayer(i, 0), h_snow_true[i], m_I0);
     Q_sw_ia[i] = -sw_in[i] * (1. - albedoValue) * (1. - i0);
-    penSW[i] = sw_in[i] * (1. - albedoValue) * i0;
+    const double extinction = 0.; // TODO: Replace with de Beer's law or a module
+    penSW[i] = sw_in[i] * (1. - albedoValue) * i0 * (1. - extinction);
+    Q_sw_base[i] = sw_in[i] * (1. - albedoValue) * i0 * extinction;
+
     // Longwave flux
     Q_lw_ia[i] = stefanBoltzmannLaw(tice.zIndexAndLayer(i, 0)) - lw_in[i];
     double dQlw_dT
         = 4 / kelvin(tice.zIndexAndLayer(i, 0)) * stefanBoltzmannLaw(tice.zIndexAndLayer(i, 0));
+
     // Total flux
     qia[i] = Q_lh_ia[i] + Q_sh_ia[i] + Q_sw_ia[i] + Q_lw_ia[i];
+
     // Total temperature dependence of flux
     dqia_dt[i] = dQlh_dT + dQsh_dT + dQlw_dT;
 }
diff --git a/physics/src/modules/FluxCalculationModule/include/FiniteElementFluxes.hpp b/physics/src/modules/FluxCalculationModule/include/FiniteElementFluxes.hpp
index 439e773f7..377a81312 100644
--- a/physics/src/modules/FluxCalculationModule/include/FiniteElementFluxes.hpp
+++ b/physics/src/modules/FluxCalculationModule/include/FiniteElementFluxes.hpp
@@ -1,18 +1,18 @@
 /*!
  * @file FiniteElementFluxes.hpp
  *
- * @date Apr 29, 2022
+ * @date 02 Sep 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
 #ifndef FINITEELEMENTFLUXES_HPP
 #define FINITEELEMENTFLUXES_HPP
 
-#include "include/ModelArrayRef.hpp"
 #include "include/Configured.hpp"
 #include "include/IFluxCalculation.hpp"
 #include "include/IIceAlbedo.hpp"
 #include "include/IIceOceanHeatFlux.hpp"
+#include "include/ModelArrayRef.hpp"
 #include "include/ModelComponent.hpp"
 
 namespace Nextsim {
@@ -25,7 +25,6 @@ class FiniteElementFluxes : public IFluxCalculation, public Configured<FiniteEle
         , evap(ModelArray::Type::H)
         , Q_lh_ow(ModelArray::Type::H)
         , Q_sh_ow(ModelArray::Type::H)
-        , Q_sw_ow(ModelArray::Type::H)
         , Q_lw_ow(ModelArray::Type::H)
         , Q_lh_ia(ModelArray::Type::H)
         , Q_sh_ia(ModelArray::Type::H)
@@ -89,7 +88,6 @@ class FiniteElementFluxes : public IFluxCalculation, public Configured<FiniteEle
     HField evap; // Open water evaporative mass flux [kg  m⁻²]
     HField Q_lh_ow; // Open water latent heat flux [W m⁻²]
     HField Q_sh_ow; // Open water sensible heat flux [W m⁻²]
-    HField Q_sw_ow; // Open water incident shortwave radiative flux [W m⁻²]
     HField Q_lw_ow; // Open water net longwave radiative flux [W m⁻²]
     HField Q_lh_ia; // Ice latent heat flux [W m⁻²]
     HField Q_sh_ia; // Ice sensible heat flux [W m⁻²]
@@ -111,8 +109,7 @@ class FiniteElementFluxes : public IFluxCalculation, public Configured<FiniteEle
     ModelArrayRef<Protected::P_AIR> p_air;
     ModelArrayRef<Protected::WIND_SPEED> v_air;
     ModelArrayRef<Protected::H_SNOW> h_snow; // cell-averaged value
-    ModelArrayRef<Protected::HTRUE_SNOW>
-        h_snow_true; // cell-averaged value
+    ModelArrayRef<Protected::HTRUE_SNOW> h_snow_true; // cell-averaged value
     ModelArrayRef<Protected::C_ICE> cice;
     ModelArrayRef<Protected::T_ICE> tice;
     ModelArrayRef<Protected::SW_IN> sw_in;
diff --git a/physics/src/modules/OceanBoundaryModule/ConfiguredOcean.cpp b/physics/src/modules/OceanBoundaryModule/ConfiguredOcean.cpp
index 03aed3154..af41440c5 100644
--- a/physics/src/modules/OceanBoundaryModule/ConfiguredOcean.cpp
+++ b/physics/src/modules/OceanBoundaryModule/ConfiguredOcean.cpp
@@ -1,7 +1,7 @@
 /*!
  * @file ConfiguredOcean.cpp
  *
- * @date 7 Sep 2023
+ * @date 30 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -110,9 +110,11 @@ void ConfiguredOcean::updateBefore(const TimestepTime& tst)
 
 void ConfiguredOcean::updateAfter(const TimestepTime& tst)
 {
+    overElements(
+        std::bind(&IOceanBoundary::mergeFluxes, this, std::placeholders::_1, std::placeholders::_2),
+        tst);
     slabOcean.update(tst);
     sst = ModelArrayRef<Protected::SLAB_SST, RO>(getStore()).data();
     sss = ModelArrayRef<Protected::SLAB_SSS, RO>(getStore()).data();
-
 }
 } /* namespace Nextsim */
diff --git a/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp b/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp
index a38a11d0b..d955cdccd 100644
--- a/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp
+++ b/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp
@@ -1,14 +1,14 @@
 /*!
  * @file TOPAZOcean.cpp
  *
- * @date 7 Sep 2023
+ * @date 30 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
 #include "include/TOPAZOcean.hpp"
 
-#include "include/IIceOceanHeatFlux.hpp"
 #include "include/IFreezingPoint.hpp"
+#include "include/IIceOceanHeatFlux.hpp"
 #include "include/Module.hpp"
 #include "include/ParaGridIO.hpp"
 #include "include/constants.hpp"
@@ -49,7 +49,6 @@ void TOPAZOcean::configure()
 
     getStore().registerArray(Protected::EXT_SST, &sstExt, RO);
     getStore().registerArray(Protected::EXT_SSS, &sssExt, RO);
-
 }
 
 void TOPAZOcean::updateBefore(const TimestepTime& tst)
@@ -65,22 +64,23 @@ void TOPAZOcean::updateBefore(const TimestepTime& tst)
     v = state.data.at("v");
 
     cpml = Water::rho * Water::cp * mld;
-    overElements(std::bind(&TOPAZOcean::updateTf, this, std::placeholders::_1,
-                     std::placeholders::_2),
+    overElements(
+        std::bind(&TOPAZOcean::updateTf, this, std::placeholders::_1, std::placeholders::_2),
         TimestepTime());
 
     Module::getImplementation<IIceOceanHeatFlux>().update(tst);
-
 }
 
 void TOPAZOcean::updateAfter(const TimestepTime& tst)
 {
+    overElements(
+        std::bind(&IOceanBoundary::mergeFluxes, this, std::placeholders::_1, std::placeholders::_2),
+        tst);
     slabOcean.update(tst);
     sst = ModelArrayRef<Protected::SLAB_SST, RO>(getStore()).data();
     sss = ModelArrayRef<Protected::SLAB_SSS, RO>(getStore()).data();
 }
 
-
 void TOPAZOcean::setFilePath(const std::string& filePathIn) { filePath = filePathIn; }
 
 void TOPAZOcean::setData(const ModelState::DataMap& ms)
diff --git a/physics/src/modules/include/IAtmosphereBoundary.hpp b/physics/src/modules/include/IAtmosphereBoundary.hpp
index 5b64a734a..ed1d7280e 100644
--- a/physics/src/modules/include/IAtmosphereBoundary.hpp
+++ b/physics/src/modules/include/IAtmosphereBoundary.hpp
@@ -1,7 +1,7 @@
 /*!
  * @file IAtmosphereBoundary.hpp
  *
- * @date Sep 22, 2022
+ * @date 30 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -15,12 +15,12 @@
 namespace Nextsim {
 
 namespace CouplingFields {
-constexpr TextTag SUBL = "SUBL"; // sublimation mass flux kg s⁻¹ m⁻²
-constexpr TextTag SNOW = "SNOW"; // snowfall mass flux kg s⁻¹ m⁻²
-constexpr TextTag RAIN = "RAIN"; // rainfall mass flux kg s⁻¹ m⁻²
-constexpr TextTag EVAP = "EVAP"; // evaporation mass flux kg s⁻¹ m⁻²
-constexpr TextTag WIND_U = "WIND_U"; // x-aligned wind component m s⁻¹
-constexpr TextTag WIND_V = "WIND_V"; // y-aligned wind component m s⁻¹
+    constexpr TextTag SUBL = "SUBL"; // sublimation mass flux kg s⁻¹ m⁻²
+    constexpr TextTag SNOW = "SNOW"; // snowfall mass flux kg s⁻¹ m⁻²
+    constexpr TextTag RAIN = "RAIN"; // rainfall mass flux kg s⁻¹ m⁻²
+    constexpr TextTag EVAP = "EVAP"; // evaporation mass flux kg s⁻¹ m⁻²
+    constexpr TextTag WIND_U = "WIND_U"; // x-aligned wind component m s⁻¹
+    constexpr TextTag WIND_V = "WIND_V"; // y-aligned wind component m s⁻¹
 
 }
 //! An interface class for the atmospheric inputs into the ice physics.
@@ -29,15 +29,17 @@ class IAtmosphereBoundary : public ModelComponent {
     IAtmosphereBoundary()
         : qia(ModelArray::Type::H)
         , dqia_dt(ModelArray::Type::H)
-    , qow(ModelArray::Type::H)
-    , subl(ModelArray::Type::H)
-    , snow(ModelArray::Type::H)
-    , rain(ModelArray::Type::H)
-    , evap(ModelArray::Type::H)
-    , emp(ModelArray::Type::H)
-    , uwind(ModelArray::Type::U)
-    , vwind(ModelArray::Type::V)
-    , penSW(ModelArray::Type::H)
+        , qow(ModelArray::Type::H)
+        , subl(ModelArray::Type::H)
+        , snow(ModelArray::Type::H)
+        , rain(ModelArray::Type::H)
+        , evap(ModelArray::Type::H)
+        , emp(ModelArray::Type::H)
+        , uwind(ModelArray::Type::U)
+        , vwind(ModelArray::Type::V)
+        , penSW(ModelArray::Type::H)
+        , qswow(ModelArray::Type::H)
+        , qswBase(ModelArray::Type::H)
     {
         m_couplingArrays.registerArray(CouplingFields::SUBL, &subl, RW);
         m_couplingArrays.registerArray(CouplingFields::SNOW, &snow, RW);
@@ -55,6 +57,8 @@ class IAtmosphereBoundary : public ModelComponent {
         getStore().registerArray(Protected::WIND_U, &uwind, RO);
         getStore().registerArray(Protected::WIND_V, &vwind, RO);
         getStore().registerArray(Shared::Q_PEN_SW, &penSW, RW);
+        getStore().registerArray(Shared::Q_SW_OW, &qswow, RW);
+        getStore().registerArray(Shared::Q_SW_BASE, &qswBase, RW);
     }
     virtual ~IAtmosphereBoundary() = default;
 
@@ -75,11 +79,12 @@ class IAtmosphereBoundary : public ModelComponent {
         uwind.resize();
         vwind.resize();
         penSW.resize();
+        qswow.resize();
+        qswBase.resize();
     }
     virtual void update(const TimestepTime& tst) { }
 
 protected:
-
     ModelArrayReferenceStore& couplingArrays() { return m_couplingArrays; }
 
     HField qia;
@@ -93,6 +98,8 @@ class IAtmosphereBoundary : public ModelComponent {
     UField uwind;
     VField vwind;
     HField penSW;
+    HField qswow;
+    HField qswBase;
 
     ModelArrayReferenceStore m_couplingArrays;
 };
diff --git a/physics/src/modules/include/IFluxCalculation.hpp b/physics/src/modules/include/IFluxCalculation.hpp
index 971247db0..0e07aa383 100644
--- a/physics/src/modules/include/IFluxCalculation.hpp
+++ b/physics/src/modules/include/IFluxCalculation.hpp
@@ -1,15 +1,15 @@
 /*!
  * @file IFluxCalculation.hpp
  *
- * @date Apr 29, 2022
+ * @date 30 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
 #ifndef IFLUXCALCULATION_HPP
 #define IFLUXCALCULATION_HPP
 
-#include "include/ModelArrayRef.hpp"
 #include "include/Configured.hpp"
+#include "include/ModelArrayRef.hpp"
 #include "include/ModelComponent.hpp"
 #include "include/ModelState.hpp"
 
@@ -23,6 +23,8 @@ class IFluxCalculation : public ModelComponent {
         , qia(getStore())
         , penSW(getStore())
         , dqia_dt(getStore())
+        , Q_sw_ow(getStore())
+        , Q_sw_base(getStore())
     {
     }
     virtual ~IFluxCalculation() = default;
@@ -53,6 +55,8 @@ class IFluxCalculation : public ModelComponent {
     ModelArrayRef<Shared::Q_IA, RW> qia;
     ModelArrayRef<Shared::Q_PEN_SW, RW> penSW;
     ModelArrayRef<Shared::DQIA_DT, RW> dqia_dt;
+    ModelArrayRef<Shared::Q_SW_OW, RW> Q_sw_ow;
+    ModelArrayRef<Shared::Q_SW_BASE, RW> Q_sw_base;
 };
 }
 #endif /* IFLUXCALCULATION_HPP */
diff --git a/physics/src/modules/include/IIceThermodynamics.hpp b/physics/src/modules/include/IIceThermodynamics.hpp
index 3546acee3..073d63c2d 100644
--- a/physics/src/modules/include/IIceThermodynamics.hpp
+++ b/physics/src/modules/include/IIceThermodynamics.hpp
@@ -1,16 +1,16 @@
 /*!
  * @file IIceThermodynamics.hpp
  *
- * @date Mar 16, 2022
+ * @date 30 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
 #ifndef IICETHERMODYNAMICS_HPP
 #define IICETHERMODYNAMICS_HPP
 
-#include "include/ModelArrayRef.hpp"
 #include "include/ConfigurationHelp.hpp"
 #include "include/ModelArray.hpp"
+#include "include/ModelArrayRef.hpp"
 #include "include/ModelComponent.hpp"
 #include "include/Time.hpp"
 
@@ -64,6 +64,7 @@ class IIceThermodynamics : public ModelComponent {
         , tf(getStore())
         , snowfall(getStore())
         , sss(getStore())
+        , qswBase(getStore())
     {
         registerModule();
 
@@ -76,6 +77,7 @@ class IIceThermodynamics : public ModelComponent {
     ModelArrayRef<Shared::H_SNOW, RW> hsnow; // From Ice Growth
     ModelArrayRef<Shared::Q_IC, RW>
         qic; // From IceTemperature. Conductive heat flux to the ice surface.
+    ModelArrayRef<Shared::Q_SW_BASE, RW> qswBase; // Short-wave flux through the base of the ice
     ModelArrayRef<Shared::Q_IO, RW> qio; // From FluxCalculation
     ModelArrayRef<Shared::Q_IA, RO> qia; // From FluxCalculation
     ModelArrayRef<Shared::DQIA_DT, RO> dQia_dt; // From FluxCalculation
diff --git a/physics/src/modules/include/IOceanBoundary.hpp b/physics/src/modules/include/IOceanBoundary.hpp
index b85c4cb06..b201608ce 100644
--- a/physics/src/modules/include/IOceanBoundary.hpp
+++ b/physics/src/modules/include/IOceanBoundary.hpp
@@ -1,7 +1,7 @@
 /*!
  * @file IOceanBoundary.hpp
  *
- * @date Sep 12, 2022
+ * @date 30 Aug 2024
  * @author Tim Spain <timothy.spain@nersc.no>
  */
 
@@ -9,20 +9,29 @@
 #define IOCEANBOUNDARY_HPP
 
 #include "include/ModelComponent.hpp"
+#include "include/constants.hpp"
 
 namespace Nextsim {
 
 namespace CouplingFields {
-constexpr TextTag SST = "SST"; // sea surface temperature ˚C
-constexpr TextTag SSS = "SSS"; // sea surface salinity PSU
-constexpr TextTag MLD = "MLD"; // Mixed layer or slab ocean depth m
-constexpr TextTag OCEAN_U = "U"; // x(east)-ward ocean current m s⁻¹
-constexpr TextTag OCEAN_V = "V"; // y(north)-ward ocean current m s⁻¹
+    constexpr TextTag SST = "SST"; // sea surface temperature ˚C
+    constexpr TextTag SSS = "SSS"; // sea surface salinity PSU
+    constexpr TextTag MLD = "MLD"; // Mixed layer or slab ocean depth m
+    constexpr TextTag OCEAN_U = "U"; // x(east)-ward ocean current m s⁻¹
+    constexpr TextTag OCEAN_V = "V"; // y(north)-ward ocean current m s⁻¹
 }
 //! An interface class for the oceanic inputs into the ice physics.
 class IOceanBoundary : public ModelComponent {
 public:
     IOceanBoundary()
+        : cice(getStore())
+        , qswow(getStore())
+        , emp(getStore())
+        , newIce(getStore())
+        , deltaHice(getStore())
+        , deltaSmelt(getStore())
+        , qow(getStore())
+        , qswBase(getStore())
     {
         m_couplingArrays.registerArray(CouplingFields::SST, &sst, RW);
         m_couplingArrays.registerArray(CouplingFields::SSS, &sss, RW);
@@ -37,6 +46,8 @@ class IOceanBoundary : public ModelComponent {
         getStore().registerArray(Protected::TF, &tf, RO);
         getStore().registerArray(Protected::OCEAN_U, &u, RO);
         getStore().registerArray(Protected::OCEAN_V, &v, RO);
+        getStore().registerArray(Protected::FWFLUX, &fwFlux, RO);
+        getStore().registerArray(Protected::SFLUX, &sFlux, RO);
     }
     virtual ~IOceanBoundary() = default;
 
@@ -54,6 +65,10 @@ class IOceanBoundary : public ModelComponent {
         tf.resize();
         u.resize();
         v.resize();
+        qNoSun.resize();
+        qswNet.resize();
+        fwFlux.resize();
+        sFlux.resize();
 
         if (ms.count("sst")) {
             sst = ms.at("sst");
@@ -77,6 +92,32 @@ class IOceanBoundary : public ModelComponent {
      */
     virtual void updateAfter(const TimestepTime& tst) = 0;
 
+    /*!
+     * Merges the ice-ocean fluxes and ocean-atmosphere fluxes into a single field to be passed to a
+     * slab-ocean implementation or an ocean model through a coupler.
+     */
+    void mergeFluxes(size_t i, const TimestepTime& tst)
+    {
+        const double dt = tst.step.seconds();
+
+        qswNet(i) = cice(i) * qswBase(i) + (1 - cice(i)) * qswow(i);
+        qNoSun(i) = cice(i) * qio(i) + (1 - cice(i)) * qow(i) - qswNet(i);
+
+        // ice volume change, both laterally and vertically
+        const double deltaIceVol = newIce(i) + deltaHice(i) * cice(i);
+        // change in snow volume due to melting (should be < 0)
+        const double meltSnowVol = deltaSmelt(i) * cice(i);
+        // Effective ice salinity is always less than or equal to the SSS, and here we use the right
+        // units too
+        const double effectiveIceSal = 1e-3 * std::min(sss(i), Ice::s);
+
+        // Positive flux is up!
+        fwFlux(i)
+            = ((1 - effectiveIceSal) * Ice::rho * deltaIceVol + Ice::rhoSnow * meltSnowVol) / dt
+            + emp(i) * (1 - cice(i));
+        sFlux(i) = effectiveIceSal * Ice::rho * deltaIceVol / dt;
+    }
+
 protected:
     HField qio; // Ice-ocean heat flux, W m⁻²
     HField sst; // Coupled or slab ocean sea surface temperature, ˚C
@@ -86,8 +127,21 @@ class IOceanBoundary : public ModelComponent {
     HField cpml; // Heat capacity of the mixed layer, J K⁻¹ m²
     UField u; // x(east)-ward ocean current, m s⁻¹
     VField v; // y(north)-ward ocean current, m s⁻¹
+    HField qNoSun; // Net surface ocean heat flux, except short wave, W m⁻²
+    HField qswNet; // Net surface ocean shortwave flux, W m⁻²
+    HField fwFlux; // Net surface ocean fresh-water flux, kg m⁻²
+    HField sFlux; // Net surface ocean salt flux, kg m⁻²
 
     ModelArrayReferenceStore m_couplingArrays;
+
+    ModelArrayRef<Protected::C_ICE, RO> cice;
+    ModelArrayRef<Shared::Q_SW_OW, RO> qswow;
+    ModelArrayRef<Protected::EVAP_MINUS_PRECIP, RO> emp;
+    ModelArrayRef<Shared::NEW_ICE, RW> newIce;
+    ModelArrayRef<Shared::DELTA_HICE, RW> deltaHice;
+    ModelArrayRef<Shared::HSNOW_MELT, RW> deltaSmelt;
+    ModelArrayRef<Shared::Q_OW, RW> qow;
+    ModelArrayRef<Shared::Q_SW_BASE, RW> qswBase;
 };
 } /* namespace Nextsim */