From f69e10113c5b1f20105546426abdeb263a53f71a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20=C3=96rn=20=C3=93lason?= Date: Fri, 30 Aug 2024 13:43:11 +0200 Subject: [PATCH] Define fluxes needed for the ocean coupling In particular, these are: the non-solar heat flux (Q_NO_SUN), the short-wave flux (Q_SW_OCN), the fresh-water flux (FWFLUX), and the salt flux (SFLUX), all at the ocean surface. I also defined the short wave flux through the base of the ice, but it's zero for now. This is an initial commit. It compiles, but is untested and still a bit messy. --- core/src/include/ModelComponent.hpp | 7 ++- physics/src/SlabOcean.cpp | 30 +++++------ physics/src/include/SlabOcean.hpp | 12 ++--- .../FiniteElementFluxes.cpp | 10 +++- .../include/FiniteElementFluxes.hpp | 2 - .../OceanBoundaryModule/ConfiguredOcean.cpp | 3 ++ .../OceanBoundaryModule/TOPAZOcean.cpp | 4 ++ .../modules/include/IAtmosphereBoundary.hpp | 8 +++ .../src/modules/include/IFluxCalculation.hpp | 4 ++ .../modules/include/IIceThermodynamics.hpp | 2 + .../src/modules/include/IOceanBoundary.hpp | 54 +++++++++++++++++++ 11 files changed, 107 insertions(+), 29 deletions(-) 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 * @author Einar Ólason */ @@ -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 ecc5d6fd6..5f7ba7a0d 100644 --- a/physics/src/SlabOcean.cpp +++ b/physics/src/SlabOcean.cpp @@ -74,33 +74,29 @@ std::unordered_set 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 */ @@ -38,9 +38,8 @@ class SlabOcean : public ModelComponent, public Configured { , 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 { ModelArrayRef qio; ModelArrayRef qow; // TODO ModelArrayRef to assimilation flux - ModelArrayRef newIce; - ModelArrayRef deltaHice; - ModelArrayRef deltaSmelt; + ModelArrayRef fwFlux; + ModelArrayRef 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 a1a6d7ec5..a1c65945e 100644 --- a/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp +++ b/physics/src/modules/FluxCalculationModule/FiniteElementFluxes.cpp @@ -59,7 +59,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(); @@ -142,22 +141,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 e243f8775..d73b7cd23 100644 --- a/physics/src/modules/FluxCalculationModule/include/FiniteElementFluxes.hpp +++ b/physics/src/modules/FluxCalculationModule/include/FiniteElementFluxes.hpp @@ -25,7 +25,6 @@ class FiniteElementFluxes : public IFluxCalculation, public Configured(getStore()).data(); sss = ModelArrayRef(getStore()).data(); diff --git a/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp b/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp index 70eb86bae..9bac61f49 100644 --- a/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp +++ b/physics/src/modules/OceanBoundaryModule/TOPAZOcean.cpp @@ -10,6 +10,7 @@ #include "include/Finalizer.hpp" #include "include/IIceOceanHeatFlux.hpp" #include "include/IFreezingPoint.hpp" +#include "include/IIceOceanHeatFlux.hpp" #include "include/NextsimModule.hpp" #include "include/ParaGridIO.hpp" #include "include/constants.hpp" @@ -76,6 +77,9 @@ void TOPAZOcean::updateBefore(const TimestepTime& 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(getStore()).data(); sss = ModelArrayRef(getStore()).data(); diff --git a/physics/src/modules/include/IAtmosphereBoundary.hpp b/physics/src/modules/include/IAtmosphereBoundary.hpp index cf4bb2731..f1cd2cf14 100644 --- a/physics/src/modules/include/IAtmosphereBoundary.hpp +++ b/physics/src/modules/include/IAtmosphereBoundary.hpp @@ -38,6 +38,8 @@ class IAtmosphereBoundary : public ModelComponent { , 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,6 +79,8 @@ class IAtmosphereBoundary : public ModelComponent { uwind.resize(); vwind.resize(); penSW.resize(); + qswow.resize(); + qswBase.resize(); } virtual void update(const TimestepTime& tst) { } @@ -92,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 4e58bd6a7..9573ba9a8 100644 --- a/physics/src/modules/include/IFluxCalculation.hpp +++ b/physics/src/modules/include/IFluxCalculation.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 qia; ModelArrayRef penSW; ModelArrayRef dqia_dt; + ModelArrayRef Q_sw_ow; + ModelArrayRef Q_sw_base; }; } #endif /* IFLUXCALCULATION_HPP */ diff --git a/physics/src/modules/include/IIceThermodynamics.hpp b/physics/src/modules/include/IIceThermodynamics.hpp index c819a5f84..f0bcd88f4 100644 --- a/physics/src/modules/include/IIceThermodynamics.hpp +++ b/physics/src/modules/include/IIceThermodynamics.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 hsnow; // From Ice Growth ModelArrayRef qic; // From IceTemperature. Conductive heat flux to the ice surface. + ModelArrayRef qswBase; // Short-wave flux through the base of the ice ModelArrayRef qio; // From FluxCalculation ModelArrayRef qia; // From FluxCalculation ModelArrayRef dQia_dt; // From FluxCalculation diff --git a/physics/src/modules/include/IOceanBoundary.hpp b/physics/src/modules/include/IOceanBoundary.hpp index e37336b8c..7cd4b2b51 100644 --- a/physics/src/modules/include/IOceanBoundary.hpp +++ b/physics/src/modules/include/IOceanBoundary.hpp @@ -9,6 +9,7 @@ #define IOCEANBOUNDARY_HPP #include "include/ModelComponent.hpp" +#include "include/constants.hpp" namespace Nextsim { @@ -23,6 +24,14 @@ namespace CouplingFields { 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 cice; + ModelArrayRef qswow; + ModelArrayRef emp; + ModelArrayRef newIce; + ModelArrayRef deltaHice; + ModelArrayRef deltaSmelt; + ModelArrayRef qow; + ModelArrayRef qswBase; }; } /* namespace Nextsim */