diff --git a/amr-wind/CFDSim.H b/amr-wind/CFDSim.H
index 12a7299e95..a69da11deb 100644
--- a/amr-wind/CFDSim.H
+++ b/amr-wind/CFDSim.H
@@ -20,6 +20,10 @@ class OversetManager;
 class ExtSolverMgr;
 class HelicsStorage;
 
+namespace transport {
+class TransportModel;
+}
+
 namespace turbulence {
 class TurbulenceModel;
 }
@@ -74,6 +78,12 @@ public:
         return m_physics_mgr.objects();
     }
 
+    transport::TransportModel& transport_model() { return *m_transport; }
+    const transport::TransportModel& transport_model() const
+    {
+        return *m_transport;
+    }
+
     turbulence::TurbulenceModel& turbulence_model() { return *m_turbulence; }
     const turbulence::TurbulenceModel& turbulence_model() const
     {
@@ -103,6 +113,9 @@ public:
 
     bool has_overset() const;
 
+    //! Instantiate the transport model based on user inputs
+    void create_transport_model();
+
     //! Instantiate the turbulence model based on user inputs
     void create_turbulence_model();
 
@@ -130,6 +143,8 @@ private:
 
     PhysicsMgr m_physics_mgr;
 
+    std::unique_ptr<transport::TransportModel> m_transport;
+
     std::unique_ptr<turbulence::TurbulenceModel> m_turbulence;
 
     std::unique_ptr<IOManager> m_io_mgr;
diff --git a/amr-wind/CFDSim.cpp b/amr-wind/CFDSim.cpp
index 0ef2b113f4..1d46cbe870 100644
--- a/amr-wind/CFDSim.cpp
+++ b/amr-wind/CFDSim.cpp
@@ -1,4 +1,5 @@
 #include "amr-wind/CFDSim.H"
+#include "amr-wind/transport_models/TransportModel.H"
 #include "amr-wind/turbulence/TurbulenceModel.H"
 #include "amr-wind/utilities/IOManager.H"
 #include "amr-wind/utilities/PostProcessing.H"
@@ -22,6 +23,17 @@ CFDSim::CFDSim(amrex::AmrCore& mesh)
 
 CFDSim::~CFDSim() = default;
 
+void CFDSim::create_transport_model()
+{
+    std::string transport_model_name = "ConstTransport";
+    {
+        amrex::ParmParse pp("transport");
+        pp.query("model", transport_model_name);
+    }
+    m_transport =
+        transport::TransportModel::create(transport_model_name, *this);
+}
+
 void CFDSim::create_turbulence_model()
 {
     std::string transport_model_name = "ConstTransport";
diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H
index 7714636850..258299dd92 100644
--- a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H
+++ b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.H
@@ -47,7 +47,7 @@ private:
     amrex::Real m_ref_theta{300.0};
 
     //! Transport model
-    std::unique_ptr<transport::TransportModel> m_transport;
+    const transport::TransportModel& m_transport;
 
     //! Thermal expansion coefficient
     std::unique_ptr<ScratchField> m_beta;
diff --git a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp
index 97bc2cb351..9d5da94fa2 100644
--- a/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp
+++ b/amr-wind/equation_systems/icns/source_terms/ABLMeanBoussinesq.cpp
@@ -17,7 +17,9 @@ namespace amr_wind::pde::icns {
  *  - `read_temperature_profile`
  *  - `tprofile_filename`
  */
-ABLMeanBoussinesq::ABLMeanBoussinesq(const CFDSim& sim) : m_mesh(sim.mesh())
+ABLMeanBoussinesq::ABLMeanBoussinesq(const CFDSim& sim)
+    : m_mesh(sim.mesh()), m_transport(sim.transport_model())
+
 {
 
     const auto& abl = sim.physics_manager().get<amr_wind::ABL>();
@@ -26,13 +28,7 @@ ABLMeanBoussinesq::ABLMeanBoussinesq(const CFDSim& sim) : m_mesh(sim.mesh())
     amrex::ParmParse pp_boussinesq_buoyancy("BoussinesqBuoyancy");
     pp_boussinesq_buoyancy.get("reference_temperature", m_ref_theta);
 
-    std::string transport_model_name = "ConstTransport";
-    {
-        amrex::ParmParse pp("transport");
-        pp.query("model", transport_model_name);
-    }
-    m_transport = transport::TransportModel::create(transport_model_name, sim);
-    m_beta = m_transport->beta();
+    m_beta = m_transport.beta();
 
     // gravity in `incflo` namespace
     amrex::ParmParse pp_incflo("incflo");
diff --git a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H
index d8153f5f13..205a24f478 100644
--- a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H
+++ b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.H
@@ -34,8 +34,9 @@ private:
     const Field& m_temperature;
 
     amrex::Vector<amrex::Real> m_gravity{0.0, 0.0, -9.81};
+
     //! Transport model
-    std::unique_ptr<transport::TransportModel> m_transport;
+    const transport::TransportModel& m_transport;
 
     //! Thermal expansion coefficient
     std::unique_ptr<ScratchField> m_beta;
diff --git a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp
index ffe4c42d08..111e2b1a35 100644
--- a/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp
+++ b/amr-wind/equation_systems/icns/source_terms/BoussinesqBuoyancy.cpp
@@ -16,14 +16,9 @@ namespace amr_wind::pde::icns {
  */
 BoussinesqBuoyancy::BoussinesqBuoyancy(const CFDSim& sim)
     : m_temperature(sim.repo().get_field("temperature"))
+    , m_transport(sim.transport_model())
 {
-    std::string transport_model_name = "ConstTransport";
-    {
-        amrex::ParmParse pp("transport");
-        pp.query("model", transport_model_name);
-    }
-    m_transport = transport::TransportModel::create(transport_model_name, sim);
-    m_beta = m_transport->beta();
+    m_beta = m_transport.beta();
 
     // gravity in `incflo` namespace
     amrex::ParmParse pp_incflo("incflo");
@@ -39,7 +34,7 @@ void BoussinesqBuoyancy::operator()(
     const FieldState fstate,
     const amrex::Array4<amrex::Real>& src_term) const
 {
-    const amrex::Real T0 = m_transport->reference_temperature();
+    const amrex::Real T0 = m_transport.reference_temperature();
     const amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> gravity{
         m_gravity[0], m_gravity[1], m_gravity[2]};
 
diff --git a/amr-wind/incflo.cpp b/amr-wind/incflo.cpp
index 9139abc7a1..64b8c53e0e 100644
--- a/amr-wind/incflo.cpp
+++ b/amr-wind/incflo.cpp
@@ -425,6 +425,7 @@ void incflo::init_physics_and_pde()
         amrex::ParmParse pp("incflo");
         pp.query("prescribe_velocity", m_prescribe_vel);
     }
+    m_sim.create_transport_model();
     m_sim.create_turbulence_model();
 
     // Initialize the refinement criteria
diff --git a/amr-wind/transport_models/ConstTransport.H b/amr-wind/transport_models/ConstTransport.H
index e64afe36cc..0117128700 100644
--- a/amr-wind/transport_models/ConstTransport.H
+++ b/amr-wind/transport_models/ConstTransport.H
@@ -152,7 +152,7 @@ public:
     }
 
     //! Return the thermal expansion coefficient
-    inline std::unique_ptr<ScratchField> beta() override
+    inline std::unique_ptr<ScratchField> beta() const override
     {
         auto beta = m_repo.create_scratch_field(1, 1);
 
@@ -184,7 +184,7 @@ public:
         return beta;
     }
 
-    inline amrex::Real reference_temperature() override
+    inline amrex::Real reference_temperature() const override
     {
         return m_reference_temperature;
     }
diff --git a/amr-wind/transport_models/TransportModel.H b/amr-wind/transport_models/TransportModel.H
index de96ba995c..f62ebdd3d3 100644
--- a/amr-wind/transport_models/TransportModel.H
+++ b/amr-wind/transport_models/TransportModel.H
@@ -44,10 +44,10 @@ public:
     scalar_diffusivity(const std::string& scalar_name) = 0;
 
     //! Thermal expansion coefficient
-    virtual std::unique_ptr<ScratchField> beta() = 0;
+    virtual std::unique_ptr<ScratchField> beta() const = 0;
 
     //! Reference temperature
-    virtual amrex::Real reference_temperature() = 0;
+    virtual amrex::Real reference_temperature() const = 0;
 };
 } // namespace amr_wind::transport
 
diff --git a/amr-wind/transport_models/TwoPhaseTransport.H b/amr-wind/transport_models/TwoPhaseTransport.H
index aec76ff854..624b0041d8 100644
--- a/amr-wind/transport_models/TwoPhaseTransport.H
+++ b/amr-wind/transport_models/TwoPhaseTransport.H
@@ -286,7 +286,7 @@ public:
     }
 
     //! Return the thermal expansion coefficient
-    inline std::unique_ptr<ScratchField> beta() override
+    inline std::unique_ptr<ScratchField> beta() const override
     {
         auto beta = m_repo.create_scratch_field(1, 1);
 
@@ -341,7 +341,7 @@ public:
         return beta;
     }
 
-    inline amrex::Real reference_temperature() override
+    inline amrex::Real reference_temperature() const override
     {
         return m_reference_temperature;
     }
diff --git a/unit_tests/offshore_wind/test_abloffshore_src.cpp b/unit_tests/offshore_wind/test_abloffshore_src.cpp
index 2fd4a91bf8..5f6eab085f 100644
--- a/unit_tests/offshore_wind/test_abloffshore_src.cpp
+++ b/unit_tests/offshore_wind/test_abloffshore_src.cpp
@@ -305,6 +305,7 @@ TEST_F(ABLOffshoreMeshTest, boussinesq)
     pde_mgr.register_icns();
     pde_mgr.register_transport_pde("Temperature");
     sim().init_physics();
+    sim().create_transport_model();
     auto& mphase = sim().physics_manager().get<amr_wind::MultiPhase>();
     // Make sure to read water level
     mphase.post_init_actions();
diff --git a/unit_tests/wind_energy/test_abl_src.cpp b/unit_tests/wind_energy/test_abl_src.cpp
index 08e4568609..4cd5698fb3 100644
--- a/unit_tests/wind_energy/test_abl_src.cpp
+++ b/unit_tests/wind_energy/test_abl_src.cpp
@@ -536,6 +536,7 @@ TEST_F(ABLMeshTest, boussinesq)
     pde_mgr.register_icns();
     pde_mgr.register_transport_pde("Temperature");
     sim().init_physics();
+    sim().create_transport_model();
 
     amr_wind::pde::icns::BoussinesqBuoyancy bb(sim());
 
@@ -582,6 +583,7 @@ TEST_F(ABLMeshTest, boussinesq_nph)
     pde_mgr.register_icns();
     pde_mgr.register_transport_pde("Temperature");
     sim().init_physics();
+    sim().create_transport_model();
 
     amr_wind::pde::icns::BoussinesqBuoyancy bb(sim());