diff --git a/CHANGELOG.md b/CHANGELOG.md
index 995fdd707b1d4547e91f0fb8b35f8ba80a9c4e7f..8b727bebfb4160afa21515317ad19706fb4a29f0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -69,15 +69,12 @@ template<class TypeTag>
 struct EnableThermalDispersion<TypeTag, TTag::MyTest> { static constexpr bool value = true; };
 ```
 
-To determine the type of dispersion tensor, please define the property `DispersionTensorType` within your `properties.hh` header. For example:
+To determine the type of dispersion tensor, please define the property `CompositionalDispersionModel` within your `properties.hh` header. Per default, the `ThermalDispersionModel` is set to the same model as the `CompositionalDispersionModel`, but this can be set specifically as well. For example:
 ```cpp
 template<class TypeTag>
-struct DispersionTensorType<TypeTag, TTag::MyTest> { using type = ScheideggersDispersionTensor<TypeTag>; };
-```
-or
-```cpp
+struct CompositionalDispersionModel<TypeTag, TTag::MyTest> { using type = ScheideggersDispersionTensor<TypeTag>; };
 template<class TypeTag>
-struct DispersionTensorType<TypeTag, TTag::MyTest> { using type = FullDispersionTensor<TypeTag>; };
+struct ThermalDispersionModel<TypeTag, TTag::MyTest> { using type = FullDispersionTensor<TypeTag>; };
 ```
 
 The parameters describing your dispersion tensor can then be included in your `spatialparameters.hh` file, and passed via input parameters. An example of this can be seen in the `test/porousmediumflow/1pnc/dispersion/` folder, and in the `test/porousmediumflow/tracer/constvel/` folders.
diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 95e1fdc9a62e946cd86a36c59d87e6f013b8da14..2a7d1fdbb015bff26031776f95e4870608b6ae17 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -150,7 +150,9 @@ struct SolutionDependentMolecularDiffusion { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
 struct HeatConductionType { using type = UndefinedProperty; };                  //!< The type for the calculation of the heat conduction fluxes
 template<class TypeTag, class MyTypeTag>
-struct DispersionTensorType { using type = UndefinedProperty; };                //!< The type for the calculation of the dispersion tensor
+struct CompositionalDispersionModel { using type = UndefinedProperty; };         //!< The type for the calculation of the compositional dispersion tensor
+template<class TypeTag, class MyTypeTag>
+struct ThermalDispersionModel { using type = UndefinedProperty; };              //!< The type for the calculation of the thermal dispersion tensor
 template<class TypeTag, class MyTypeTag>
 struct SolutionDependentHeatConduction { using type = UndefinedProperty; };     //!< specifies if the parameters for the heat conduction fluxes depend on the solution
 
diff --git a/dumux/flux/box/dispersionflux.hh b/dumux/flux/box/dispersionflux.hh
index 4607b3f0afb3b29eef98507faf98992e90873b08..bb5e81bee461f115820df68bc2c54800a73f3c69 100644
--- a/dumux/flux/box/dispersionflux.hh
+++ b/dumux/flux/box/dispersionflux.hh
@@ -38,14 +38,14 @@ namespace Dumux {
 
 // forward declaration
 template<class TypeTag, class DiscretizationMethod, ReferenceSystemFormulation referenceSystem>
-class OnePDispersionFluxImplementation;
+class DispersionFluxImplementation;
 
 /*!
  * \ingroup BoxFlux
  * \brief Specialization of a dispersion flux for the box method
  */
 template <class TypeTag, ReferenceSystemFormulation referenceSystem>
-class OnePDispersionFluxImplementation<TypeTag, DiscretizationMethods::Box, referenceSystem>
+class DispersionFluxImplementation<TypeTag, DiscretizationMethods::Box, referenceSystem>
 {
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
@@ -65,6 +65,7 @@ class OnePDispersionFluxImplementation<TypeTag, DiscretizationMethods::Box, refe
     using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
     using Element = typename GridView::template Codim<0>::Entity;
     using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+    using Indices = typename ModelTraits::Indices;
 
     enum { dim = GridView::dimension} ;
     enum { dimWorld = GridView::dimensionworld} ;
@@ -102,10 +103,6 @@ public:
     {
         ComponentFluxVector componentFlux(0.0);
 
-        // collect the dispersion tensor, the fluxVarsCache and the shape values
-        const DimWorldMatrix dispersionTensor = VolumeVariables::DispersionTensorType::dispersionTensor(problem, scvf, fvGeometry,
-                                                                                                       elemVolVars, elemFluxVarsCache);
-
         const auto& fluxVarsCache = elemFluxVarsCache[scvf];
         const auto& shapeValues = fluxVarsCache.shapeValues();
 
@@ -119,6 +116,12 @@ public:
 
         for (int compIdx = 0; compIdx < numComponents; compIdx++)
         {
+            // collect the dispersion tensor, the fluxVarsCache and the shape values
+            const auto& dispersionTensor =
+                ModelTraits::CompositionalDispersionModel::compositionalDispersionTensor(problem, scvf, fvGeometry,
+                                                                                         elemVolVars, elemFluxVarsCache,
+                                                                                         phaseIdx, compIdx);
+
             // the mole/mass fraction gradient
             Dune::FieldVector<Scalar, dimWorld> gradX(0.0);
             for (auto&& scv : scvs(fvGeometry))
@@ -148,8 +151,10 @@ public:
                                                 const ElementFluxVariablesCache& elemFluxVarsCache)
     {
         // collect the dispersion tensor
-        const DimWorldMatrix dispersionTensor = VolumeVariables::DispersionTensorType::dispersionTensor(problem, scvf, fvGeometry,
-                                                                                                       elemVolVars, elemFluxVarsCache);
+        const auto& dispersionTensor =
+            ModelTraits::ThermalDispersionModel::thermalDispersionTensor(problem, scvf, fvGeometry,
+                                                                         elemVolVars, elemFluxVarsCache,
+                                                                         phaseIdx);
         // compute the temperature gradient with the shape functions
         const auto& fluxVarsCache = elemFluxVarsCache[scvf];
         Dune::FieldVector<Scalar, GridView::dimensionworld> gradTemp(0.0);
diff --git a/dumux/flux/cctpfa/dispersionflux.hh b/dumux/flux/cctpfa/dispersionflux.hh
index af2dad28aa054580e7fe297b568edecde51e8e0d..a89024c9e6250421318aa1e84b3bc0090ebde27f 100644
--- a/dumux/flux/cctpfa/dispersionflux.hh
+++ b/dumux/flux/cctpfa/dispersionflux.hh
@@ -39,14 +39,14 @@ namespace Dumux {
 
 // forward declaration
 template<class TypeTag, class DiscretizationMethod, ReferenceSystemFormulation referenceSystem>
-class OnePDispersionFluxImplementation;
+class DispersionFluxImplementation;
 
 /*!
  * \ingroup CCTpfaFlux
  * \brief Specialization of a Dispersion flux for the cctpfa method
  */
 template <class TypeTag, ReferenceSystemFormulation referenceSystem>
-class OnePDispersionFluxImplementation<TypeTag, DiscretizationMethods::CCTpfa, referenceSystem>
+class DispersionFluxImplementation<TypeTag, DiscretizationMethods::CCTpfa, referenceSystem>
 {
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
@@ -66,6 +66,7 @@ class OnePDispersionFluxImplementation<TypeTag, DiscretizationMethods::CCTpfa, r
     using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
     using Element = typename GridView::template Codim<0>::Entity;
     using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
+    using Indices = typename ModelTraits::Indices;
 
     enum { dim = GridView::dimension} ;
     enum { dimWorld = GridView::dimensionworld} ;
@@ -110,16 +111,18 @@ public:
         const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
         const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
-        const DimWorldMatrix dispersionTensor = VolumeVariables::DispersionTensorType::dispersionTensor(problem, scvf, fvGeometry,
-                                                                                                        elemVolVars, elemFluxVarsCache);
-        const auto dij = computeTpfaTransmissibility(scvf, fvGeometry.scv(scvf.insideScvIdx()), dispersionTensor, insideVolVars.extrusionFactor());
-
         const auto rhoInside = massOrMolarDensity(insideVolVars, referenceSystem, phaseIdx);
         const auto rhoOutside = massOrMolarDensity(outsideVolVars, referenceSystem, phaseIdx);
         const Scalar rho = 0.5*(rhoInside + rhoOutside);
 
         for (int compIdx = 0; compIdx < numComponents; compIdx++)
         {
+            const auto& dispersionTensor =
+                ModelTraits::CompositionalDispersionModel::compositionalDispersionTensor(problem, scvf, fvGeometry,
+                                                                                         elemVolVars, elemFluxVarsCache,
+                                                                                         phaseIdx, compIdx);
+            const auto dij = computeTpfaTransmissibility(scvf, fvGeometry.scv(scvf.insideScvIdx()), dispersionTensor, insideVolVars.extrusionFactor());
+
             const auto xInside = massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx);
             const auto xOutide = massOrMoleFraction(outsideVolVars, referenceSystem, phaseIdx, compIdx);
 
@@ -128,7 +131,6 @@ public:
         return componentFlux;
     }
 
-
     /*!
      * \brief Returns the thermal dispersive flux
      *        across the given sub-control volume face.
@@ -149,8 +151,10 @@ public:
         const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
         const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
-        const DimWorldMatrix dispersionTensor = VolumeVariables::DispersionTensorType::dispersionTensor(problem, scvf, fvGeometry,
-                                                                                                        elemVolVars, elemFluxVarsCache);
+        const auto& dispersionTensor =
+            ModelTraits::ThermalDispersionModel::thermalDispersionTensor(problem, scvf, fvGeometry,
+                                                                         elemVolVars, elemFluxVarsCache,
+                                                                         phaseIdx);
         const auto dij = computeTpfaTransmissibility(scvf, fvGeometry.scv(scvf.insideScvIdx()), dispersionTensor, insideVolVars.extrusionFactor());
 
         // get the inside/outside temperatures
diff --git a/dumux/flux/dispersionflux_fwd.hh b/dumux/flux/dispersionflux_fwd.hh
index 8112a15a0bc1d238ebb186c0777b40258785c823..fd402ea0d437c8c930e4e4b4e6b8fe7f78410600 100644
--- a/dumux/flux/dispersionflux_fwd.hh
+++ b/dumux/flux/dispersionflux_fwd.hh
@@ -31,16 +31,16 @@ namespace Dumux {
 
 // declaration of primary template
 template<class TypeTag, class DiscretizationMethod, ReferenceSystemFormulation referenceSystem>
-class OnePDispersionFluxImplementation;
+class DispersionFluxImplementation;
 
 /*!
  * \ingroup Flux
  * \brief Evaluates the dispersive flux
  */
 template<class TypeTag, ReferenceSystemFormulation referenceSystem = ReferenceSystemFormulation::massAveraged>
-using OnePDiffusiveDispersionFlux = OnePDispersionFluxImplementation<TypeTag,
-                                                                     typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod,
-                                                                     referenceSystem>;
+using DiffusiveDispersionFlux = DispersionFluxImplementation<TypeTag,
+                                                             typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod,
+                                                             referenceSystem>;
 
 } // end namespace Dumux
 
diff --git a/dumux/material/fluidmatrixinteractions/dispersiontensors/fulltensor.hh b/dumux/material/fluidmatrixinteractions/dispersiontensors/fulltensor.hh
index 851133c6434f2c81417ffd8a702d9ce2ffe86dce..911ce3a371302e9dd487c71460880e81052a5fed 100644
--- a/dumux/material/fluidmatrixinteractions/dispersiontensors/fulltensor.hh
+++ b/dumux/material/fluidmatrixinteractions/dispersiontensors/fulltensor.hh
@@ -43,12 +43,23 @@ class FullDispersionTensor
 public:
 
     template <class ElementFluxVariablesCache>
-    static DimWorldMatrix dispersionTensor([[maybe_unused]] const Problem& problem,
-                                           [[maybe_unused]] const SubControlVolumeFace& scvf,
-                                           [[maybe_unused]] const FVElementGeometry& fvGeometry,
-                                           [[maybe_unused]] const ElementVolumeVariables& elemVolVars,
-                                           [[maybe_unused]] const ElementFluxVariablesCache& elemFluxVarsCache)
-    { return problem.spatialParams().dispersionTensor(scvf.center()); }
+    static DimWorldMatrix compositionalDispersionTensor(const Problem& problem,
+                                                        const SubControlVolumeFace& scvf,
+                                                        const FVElementGeometry&,
+                                                        const ElementVolumeVariables&,
+                                                        const ElementFluxVariablesCache&,
+                                                        const int phaseIdx,
+                                                        const int compIdx)
+    { return problem.spatialParams().dispersionTensor(scvf.center(), phaseIdx, compIdx); }
+
+    template <class ElementFluxVariablesCache>
+    static DimWorldMatrix thermalDispersionTensor(const Problem& problem,
+                                                  const SubControlVolumeFace& scvf,
+                                                  const FVElementGeometry&,
+                                                  const ElementVolumeVariables&,
+                                                  const ElementFluxVariablesCache&,
+                                                  const int phaseIdx)
+    { return problem.spatialParams().dispersionTensor(scvf.center(), phaseIdx); }
 
 };
 
diff --git a/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh b/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh
index d3e1ad30f8c122380ef8a1bf5ba0f553df7215f1..23cd65c4b8117d78bd95e5de1ac1d621d0518932 100644
--- a/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh
+++ b/dumux/material/fluidmatrixinteractions/dispersiontensors/scheidegger.hh
@@ -60,14 +60,53 @@ class ScheideggersDispersionTensor
 
 public:
     template <class ElementFluxVariablesCache>
-    static DimWorldMatrix dispersionTensor(const Problem& problem,
-                                           const SubControlVolumeFace& scvf,
-                                           [[maybe_unused]] const FVElementGeometry& fvGeometry,
-                                           [[maybe_unused]] const ElementVolumeVariables& elemVolVars,
-                                           [[maybe_unused]] const ElementFluxVariablesCache& elemFluxVarsCache)
+    static DimWorldMatrix compositionalDispersionTensor(const Problem& problem,
+                                                        const SubControlVolumeFace& scvf,
+                                                        const FVElementGeometry& fvGeometry,
+                                                        const ElementVolumeVariables& elemVolVars,
+                                                        const ElementFluxVariablesCache& elemFluxVarsCache,
+                                                        const int phaseIdx,
+                                                        const int compIdx)
     {
         DimWorldMatrix dispersionTensor(0.0);
 
+        // Get the velocity either from the reconstruction, or from the spatialparams
+        auto velocity = dispersionVelocity_(problem, scvf, fvGeometry, elemVolVars, elemFluxVarsCache);
+
+        // collect the dispersion alphas at this location
+        std::array<Scalar,2> dispersivity = problem.spatialParams().dispersionAlphas(scvf.center(), phaseIdx, compIdx);
+
+        return scheideggerTensor_(dispersivity, velocity);
+    }
+
+    template <class ElementFluxVariablesCache>
+    static DimWorldMatrix thermalDispersionTensor(const Problem& problem,
+                                                  const SubControlVolumeFace& scvf,
+                                                  const FVElementGeometry& fvGeometry,
+                                                  const ElementVolumeVariables& elemVolVars,
+                                                  const ElementFluxVariablesCache& elemFluxVarsCache,
+                                                  const int phaseIdx)
+    {
+        DimWorldMatrix dispersionTensor(0.0);
+
+        // Get the velocity either from the reconstruction, or from the spatialparams
+        auto velocity = dispersionVelocity_(problem, scvf, fvGeometry, elemVolVars, elemFluxVarsCache);
+
+        // collect the dispersion alphas at this location
+        std::array<Scalar,2> dispersivity = problem.spatialParams().dispersionAlphas(scvf.center(), phaseIdx); //TODO: fix this?
+
+        return scheideggerTensor_(dispersivity, velocity);
+    }
+
+private:
+
+    template <class ElementFluxVariablesCache>
+    static Dune::FieldVector<Scalar, dimWorld> dispersionVelocity_(const Problem& problem,
+                                                                   const SubControlVolumeFace& scvf,
+                                                                   [[maybe_unused]] const FVElementGeometry& fvGeometry,
+                                                                   [[maybe_unused]] const ElementVolumeVariables& elemVolVars,
+                                                                   [[maybe_unused]] const ElementFluxVariablesCache& elemFluxVarsCache)
+    {
         // Calculate Darcy's velocity
         Dune::FieldVector<Scalar, dimWorld> velocity(0.0);
         if constexpr (stationaryVelocityField)
@@ -117,31 +156,34 @@ public:
                 DUNE_THROW(Dune::NotImplemented, "\n Scheidegger Dispersion for compositional models without given constant velocity field is only implemented using the Box method.");
         }
 
-        //calculate dispersion tensor
+        return velocity;
+    }
 
-        // collect the dispersion alphas at this location
-        std::array<Scalar,2> dispersivity = problem.spatialParams().dispersionAlphas(scvf.center());
+    static DimWorldMatrix scheideggerTensor_(const std::array<Scalar,2>& dispersivity,
+                                             const Dune::FieldVector<Scalar, dimWorld>& velocity)
+    {
+        DimWorldMatrix scheideggerTensor(0.0);
 
         //matrix multiplication of the velocity at the interface: vv^T
         for (int i=0; i < dimWorld; i++)
             for (int j = 0; j < dimWorld; j++)
-                dispersionTensor[i][j] = velocity[i]*velocity[j];
+                scheideggerTensor[i][j] = velocity[i]*velocity[j];
 
         //normalize velocity product --> vv^T/||v||, [m/s]
         Scalar vNorm = velocity.two_norm();
 
-        dispersionTensor /= vNorm;
+        scheideggerTensor /= vNorm;
         if (vNorm < 1e-20)
-            dispersionTensor = 0;
+            scheideggerTensor = 0;
 
         //multiply with dispersivity difference: vv^T/||v||*(alphaL - alphaT), [m^2/s] --> alphaL = longitudinal disp., alphaT = transverse disp.
-        dispersionTensor *= (dispersivity[0] - dispersivity[1]);
+        scheideggerTensor *= (dispersivity[0] - dispersivity[1]);
 
         //add ||v||*alphaT to the main diagonal:vv^T/||v||*(alphaL - alphaT) + ||v||*alphaT, [m^2/s]
         for (int i = 0; i < dimWorld; i++)
-            dispersionTensor[i][i] += vNorm*dispersivity[1];
+            scheideggerTensor[i][i] += vNorm*dispersivity[1];
 
-        return dispersionTensor;
+        return scheideggerTensor;
     }
 
 };
diff --git a/dumux/porenetwork/1pnc/model.hh b/dumux/porenetwork/1pnc/model.hh
index fd09f0174bcae504544ad4245d5edcd73c472671..8fe5c81fecc2f0d133f3dad1b031227ef1dbf40c 100644
--- a/dumux/porenetwork/1pnc/model.hh
+++ b/dumux/porenetwork/1pnc/model.hh
@@ -158,19 +158,17 @@ private:
     static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
     using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
 
-    using DTT = GetPropType<TypeTag, Properties::DispersionTensorType>;
     using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
     using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
-    template<class BaseTraits, class DTT, class DT, class EDM>
+    template<class BaseTraits, class DT, class EDM>
     struct NCTraits : public BaseTraits
     {
-        using DispersionTensorType = DTT;
         using DiffusionType = DT;
         using EffectiveDiffusivityModel = EDM;
     };
 
 public:
-    using type = PoreNetwork::OnePNCVolumeVariables<NCTraits<BaseTraits, DTT, DT, EDM>>;
+    using type = PoreNetwork::OnePNCVolumeVariables<NCTraits<BaseTraits, DT, EDM>>;
 };
 
 //!< Set the vtk output fields specific to this model
@@ -191,12 +189,15 @@ struct ModelTraits<TypeTag, TTag::PNMOnePNCNI>
 {
 private:
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using CompositionalDispersionModel = GetPropType<TypeTag, Properties::CompositionalDispersionModel>;
+    using ThermalDispersionModel = GetPropType<TypeTag, Properties::CompositionalDispersionModel>;
     using IsothermalTraits = OnePNCModelTraits<FluidSystem::numComponents, getPropValue<TypeTag, Properties::UseMoles>(),
                                                                            getPropValue<TypeTag, Properties::EnableCompositionalDispersion>(),
                                                                            getPropValue<TypeTag, Properties::EnableThermalDispersion>(),
-                                                                           getPropValue<TypeTag, Properties::ReplaceCompEqIdx>()>;
+                                                                           getPropValue<TypeTag, Properties::ReplaceCompEqIdx>(),
+                                                                           CompositionalDispersionModel>;
 public:
-    using type = PorousMediumFlowNIModelTraits<IsothermalTraits>;
+    using type = PorousMediumFlowNIModelTraits<IsothermalTraits, ThermalDispersionModel>;
 };
 
 template<class TypeTag>
@@ -216,21 +217,19 @@ private:
     static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
     using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
 
-    using DTT = GetPropType<TypeTag, Properties::DispersionTensorType>;
     using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
     using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
     using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
-    template<class BaseTraits, class DTT, class DT, class EDM, class ETCM>
+    template<class BaseTraits, class DT, class EDM, class ETCM>
     struct NCNITraits : public BaseTraits
     {
-        using DispersionTensorType = DTT;
         using DiffusionType = DT;
         using EffectiveDiffusivityModel = EDM;
         using EffectiveThermalConductivityModel = ETCM;
     };
 
 public:
-    using type = PoreNetwork::OnePNCVolumeVariables<NCNITraits<BaseTraits, DTT, DT, EDM, ETCM>>;
+    using type = PoreNetwork::OnePNCVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>>;
 };
 
 template<class TypeTag>
diff --git a/dumux/porousmediumflow/1pnc/model.hh b/dumux/porousmediumflow/1pnc/model.hh
index 239424c87a8f58a4bdf29f5b7ed9f5e7aad3109c..56d03bf59d3c9b009eb76d3bd76da255a1cc57fa 100644
--- a/dumux/porousmediumflow/1pnc/model.hh
+++ b/dumux/porousmediumflow/1pnc/model.hh
@@ -85,10 +85,11 @@ namespace Dumux {
  *
  * \tparam nComp the number of components to be considered.
  */
-template<int nComp, bool useM, int enableCompDisp, int enableThermDisp, int repCompEqIdx = nComp>
+template<int nComp, bool useM, int enableCompDisp, int enableThermDisp, int repCompEqIdx = nComp, class CDM = void>
 struct OnePNCModelTraits
 {
     using Indices = OnePNCIndices;
+    using CompositionalDispersionModel = CDM;
 
     static constexpr int numEq() { return nComp; }
     static constexpr int numFluidPhases() { return 1; }
@@ -130,13 +131,14 @@ struct BaseModelTraits<TypeTag, TTag::OnePNC>
 {
 private:
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-
+    using CDM = GetPropType<TypeTag, Properties::CompositionalDispersionModel>;
 public:
     using type = OnePNCModelTraits<FluidSystem::numComponents,
                                    getPropValue<TypeTag, Properties::UseMoles>(),
                                    getPropValue<TypeTag, Properties::EnableCompositionalDispersion>(),
                                    getPropValue<TypeTag, Properties::EnableThermalDispersion>(),
-                                   getPropValue<TypeTag, Properties::ReplaceCompEqIdx>()>;
+                                   getPropValue<TypeTag, Properties::ReplaceCompEqIdx>(),
+                                   CDM>;
 };
 template<class TypeTag>
 struct ModelTraits<TypeTag, TTag::OnePNC> { using type = GetPropType<TypeTag, Properties::BaseModelTraits>; }; //!< default the actually used traits to the base traits
@@ -190,19 +192,17 @@ private:
     static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
     using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
 
-    using DTT = GetPropType<TypeTag, Properties::DispersionTensorType>;
     using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
     using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
-    template<class BaseTraits, class DTT, class DT, class EDM>
+    template<class BaseTraits, class DT, class EDM>
     struct NCTraits : public BaseTraits
     {
-        using DispersionTensorType = DTT;
         using DiffusionType = DT;
         using EffectiveDiffusivityModel = EDM;
     };
 
 public:
-    using type = OnePNCVolumeVariables<NCTraits<BaseTraits, DTT, DT, EDM>>;
+    using type = OnePNCVolumeVariables<NCTraits<BaseTraits, DT, EDM>>;
 };
 
 //! Set the vtk output fields specific to this model
@@ -228,8 +228,9 @@ struct ModelTraits<TypeTag, TTag::OnePNCNI>
 {
 private:
     using IsothermalTraits = GetPropType<TypeTag, Properties::BaseModelTraits>;
+    using TDM = GetPropType<TypeTag, Properties::ThermalDispersionModel>;
 public:
-    using type = PorousMediumFlowNIModelTraits<IsothermalTraits>;
+    using type = PorousMediumFlowNIModelTraits<IsothermalTraits, TDM>;
 };
 
 template<class TypeTag>
@@ -249,21 +250,19 @@ private:
     static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
     using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
 
-    using DTT = GetPropType<TypeTag, Properties::DispersionTensorType>;
     using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
     using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
     using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
-    template<class BaseTraits, class DTT, class DT, class EDM, class ETCM>
+    template<class BaseTraits, class DT, class EDM, class ETCM>
     struct NCNITraits : public BaseTraits
     {
-        using DispersionTensorType = DTT;
         using DiffusionType = DT;
         using EffectiveDiffusivityModel = EDM;
         using EffectiveThermalConductivityModel = ETCM;
     };
 
 public:
-    using type = OnePNCVolumeVariables<NCNITraits<BaseTraits, DTT, DT, EDM, ETCM>>;
+    using type = OnePNCVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>>;
 };
 
 } // end namespace Properties
@@ -354,22 +353,20 @@ private:
     using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
     using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
 
-    using DTT = GetPropType<TypeTag, Properties::DispersionTensorType>;
     using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
     using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
     using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
-    template<class BaseTraits, class DTT, class DT, class EDM, class ETCM>
+    template<class BaseTraits, class DT, class EDM, class ETCM>
     struct NCNITraits : public BaseTraits
     {
-        using DispersionTensorType = DTT;
         using DiffusionType = DT;
         using EffectiveDiffusivityModel = EDM;
         using EffectiveThermalConductivityModel = ETCM;
     };
 
-    using EquilibriumVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, DTT, DT, EDM, ETCM>>;
+    using EquilibriumVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>>;
 public:
-    using type = NonEquilibriumVolumeVariables<NCNITraits<BaseTraits, DTT, DT, EDM, ETCM>, EquilibriumVolVars>;
+    using type = NonEquilibriumVolumeVariables<NCNITraits<BaseTraits, DT, EDM, ETCM>, EquilibriumVolVars>;
 };
 
 } // end namespace Properties
diff --git a/dumux/porousmediumflow/1pnc/volumevariables.hh b/dumux/porousmediumflow/1pnc/volumevariables.hh
index ac3a4cbd7eb7ff4628151694c3ceab9b9d39b00b..58e84c9ba00799d972111150e765a3047eb09c51 100644
--- a/dumux/porousmediumflow/1pnc/volumevariables.hh
+++ b/dumux/porousmediumflow/1pnc/volumevariables.hh
@@ -64,8 +64,6 @@ class OnePNCVolumeVariables
     };
 
 public:
-    //! Export the dispersion tensor law
-    using DispersionTensorType = typename Traits::DispersionTensorType;
     //! Export fluid state type
     using FluidState = typename Traits::FluidState;
     //! Export fluid system type
diff --git a/dumux/porousmediumflow/1pncmin/model.hh b/dumux/porousmediumflow/1pncmin/model.hh
index 21c066bbf2797d158154b559aa583f218aac061a..0db57078fc8724324726532656b1ad46f24f7d79 100644
--- a/dumux/porousmediumflow/1pncmin/model.hh
+++ b/dumux/porousmediumflow/1pncmin/model.hh
@@ -115,20 +115,20 @@ private:
     static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
     using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
 
-    using DTT = GetPropType<TypeTag, Properties::DispersionTensorType>;
+    using CDM = GetPropType<TypeTag, Properties::CompositionalDispersionModel>;
     using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
     using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
-    template<class BaseTraits, class DTT, class DT, class EDM>
+    template<class BaseTraits, class CDM, class DT, class EDM>
     struct NCTraits : public BaseTraits
     {
-        using DispersionTensorType = DTT;
+        using CompositionalDispersionModel = CDM;
         using DiffusionType = DT;
         using EffectiveDiffusivityModel = EDM;
     };
 
-    using NonMinVolVars = OnePNCVolumeVariables<NCTraits<BaseTraits, DTT, DT, EDM>>;
+    using NonMinVolVars = OnePNCVolumeVariables<NCTraits<BaseTraits, CDM, DT, EDM>>;
 public:
-    using type = MineralizationVolumeVariables<NCTraits<BaseTraits, DTT, DT, EDM>, NonMinVolVars>;
+    using type = MineralizationVolumeVariables<NCTraits<BaseTraits, CDM, DT, EDM>, NonMinVolVars>;
 };
 
 // Use the mineralization local residual
@@ -203,21 +203,21 @@ private:
     static_assert(FST::numPhases == MT::numFluidPhases(), "Number of phases mismatch between model and fluid state");
     using BaseTraits = OnePVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT>;
 
-    using DTT = GetPropType<TypeTag, Properties::DispersionTensorType>;
+    using CDM = GetPropType<TypeTag, Properties::CompositionalDispersionModel>;
     using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
     using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
     using ETCM = GetPropType< TypeTag, Properties:: ThermalConductivityModel>;
-    template<class BaseTraits, class DTT, class DT, class EDM, class ETCM>
+    template<class BaseTraits, class CDM, class DT, class EDM, class ETCM>
     struct NCNITraits : public BaseTraits
     {
-        using DispersionTensorType = DTT;
+        using CompositionalDispersionModel = CDM;
         using DiffusionType = DT;
         using EffectiveDiffusivityModel = EDM;
         using EffectiveThermalConductivityModel = ETCM;
     };
-    using NonMinVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, DTT, DT, EDM, ETCM>>;
+    using NonMinVolVars = OnePNCVolumeVariables<NCNITraits<BaseTraits, CDM, DT, EDM, ETCM>>;
 public:
-    using type = MineralizationVolumeVariables<NCNITraits<BaseTraits, DTT, DT, EDM, ETCM>, NonMinVolVars>;
+    using type = MineralizationVolumeVariables<NCNITraits<BaseTraits, CDM, DT, EDM, ETCM>, NonMinVolVars>;
 };
 //! Use the average for effective conductivities
 template<class TypeTag>
diff --git a/dumux/porousmediumflow/fluxvariables.hh b/dumux/porousmediumflow/fluxvariables.hh
index d5e40cd33de018b581aa883ef8bdeadf7ca50289..77dddf9dc2de48d6e68b11f596cf6cf51e6fcf4f 100644
--- a/dumux/porousmediumflow/fluxvariables.hh
+++ b/dumux/porousmediumflow/fluxvariables.hh
@@ -132,6 +132,7 @@ public:
     Dune::FieldVector<Scalar, numComponents> compositionalDispersionFlux([[maybe_unused]] const int phaseIdx) const
     {
         if constexpr (enableCompositionalDispersion)
+        {
             return DispersionFluxType::compositionalDispersionFlux(this->problem(),
                                                                    this->element(),
                                                                    this->fvGeometry(),
@@ -139,6 +140,7 @@ public:
                                                                    this->scvFace(),
                                                                    phaseIdx,
                                                                    this->elemFluxVarsCache());
+        }
         else
             return Dune::FieldVector<Scalar, numComponents>(0.0);
     }
@@ -149,6 +151,7 @@ public:
     Dune::FieldVector<Scalar, 1> thermalDispersionFlux([[maybe_unused]] const int phaseIdx = 0) const
     {
         if constexpr (enableThermalDispersion)
+        {
             return DispersionFluxType::thermalDispersionFlux(this->problem(),
                                                              this->element(),
                                                              this->fvGeometry(),
@@ -156,6 +159,7 @@ public:
                                                              this->scvFace(),
                                                              phaseIdx,
                                                              this->elemFluxVarsCache());
+        }
         else
             return Dune::FieldVector<Scalar, 1>(0.0);
     }
diff --git a/dumux/porousmediumflow/nonisothermal/model.hh b/dumux/porousmediumflow/nonisothermal/model.hh
index c2e5bcb6efe0054163f08c88d8eb697db02e7d36..54a76de78a31ef02d507ea7a23508e8a9d823a1f 100644
--- a/dumux/porousmediumflow/nonisothermal/model.hh
+++ b/dumux/porousmediumflow/nonisothermal/model.hh
@@ -65,12 +65,15 @@ namespace Dumux {
  *
  * \tparam IsothermalTraits Model traits of the isothermal model
  */
-template<class IsothermalT>
+template<class IsothermalT, class TDM = void>
 struct PorousMediumFlowNIModelTraits : public IsothermalT
 {
     //! Export the isothermal model traits
     using IsothermalTraits = IsothermalT;
 
+    //! Export the thermal dispersion tensor type
+    using ThermalDispersionModel = TDM;
+
     //! We solve for one more equation, i.e. the energy balance
     static constexpr int numEq() { return IsothermalTraits::numEq()+1; }
     //! only one energy equation is needed when assuming thermal equilibrium
diff --git a/dumux/porousmediumflow/properties.hh b/dumux/porousmediumflow/properties.hh
index a357494bbeb77fc33142fb16fb13f4fe5d158e0d..c362bc318fe5101414f468912cba7bf0d8fe9d3c 100644
--- a/dumux/porousmediumflow/properties.hh
+++ b/dumux/porousmediumflow/properties.hh
@@ -85,11 +85,16 @@ struct EnableThermalDispersion<TypeTag, TTag::PorousMediumFlow> { static constex
 
 //! By default, we use a diffusive flux for the dispersive fluxes
 template<class TypeTag>
-struct DispersionFluxType<TypeTag, TTag::PorousMediumFlow> { using type = OnePDiffusiveDispersionFlux<TypeTag>; };
+struct DispersionFluxType<TypeTag, TTag::PorousMediumFlow> { using type = DiffusiveDispersionFlux<TypeTag>; };
 
 //! By default, we use Scheideggers's law for the dispersive tensor calculation
 template<class TypeTag>
-struct DispersionTensorType<TypeTag, TTag::PorousMediumFlow> { using type = ScheideggersDispersionTensor<TypeTag>; };
+struct CompositionalDispersionModel<TypeTag, TTag::PorousMediumFlow> { using type = ScheideggersDispersionTensor<TypeTag>; };
+
+//! By default, we use the same dispersion tensor type as the componsitonal dispersion for the thermal disperion tensor
+template<class TypeTag>
+struct ThermalDispersionModel<TypeTag, TTag::PorousMediumFlow>
+{ using type = GetPropType<TypeTag, Properties::CompositionalDispersionModel>; };
 
 //! By default, we use fourier's law as the default for heat conduction fluxes
 template<class TypeTag>
diff --git a/dumux/porousmediumflow/tracer/model.hh b/dumux/porousmediumflow/tracer/model.hh
index 2724528e9df61d10da4eedd8fba7f6103fc5a08e..e4c7f31237e10ad72599b461c2d84ee3467b8245 100644
--- a/dumux/porousmediumflow/tracer/model.hh
+++ b/dumux/porousmediumflow/tracer/model.hh
@@ -69,10 +69,11 @@ namespace Dumux {
  * \tparam nComp the number of components to be considered.
  * \tparam useMol whether mole or mass balances are used
  */
-template<int nComp, bool useMol, bool enableCompDisp>
+template<int nComp, bool useMol, bool enableCompDisp, class CDM = void>
 struct TracerModelTraits
 {
     using Indices = TracerIndices;
+    using CompositionalDispersionModel = CDM;
 
     static constexpr int numEq() { return nComp; }
     static constexpr int numFluidPhases() { return 1; }
@@ -94,7 +95,7 @@ struct TracerModelTraits
  * \tparam FSY The fluid system type
  * \tparam MT The model traits
  */
-template<class PV, class FSY, class SSY, class SST, class MT, class DTT, class DT, class EDM>
+template<class PV, class FSY, class SSY, class SST, class MT, class DT, class EDM>
 struct TracerVolumeVariablesTraits
 {
     using PrimaryVariables = PV;
@@ -102,7 +103,6 @@ struct TracerVolumeVariablesTraits
     using SolidSystem = SSY;
     using SolidState = SST;
     using ModelTraits = MT;
-    using DispersionTensorType = DTT;
     using DiffusionType = DT;
     using EffectiveDiffusivityModel = EDM;
 };
@@ -134,10 +134,12 @@ struct ModelTraits<TypeTag, TTag::Tracer>
 {
 private:
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using CDM = GetPropType<TypeTag, Properties::CompositionalDispersionModel>;
 public:
     using type = TracerModelTraits<FluidSystem::numComponents,
                                    getPropValue<TypeTag, Properties::UseMoles>(),
-                                   getPropValue<TypeTag, Properties::EnableCompositionalDispersion>()>;
+                                   getPropValue<TypeTag, Properties::EnableCompositionalDispersion>(),
+                                   CDM>;
 };
 
 //! Use the tracer local residual function for the tracer model
@@ -158,11 +160,10 @@ private:
     using SSY = GetPropType<TypeTag, Properties::SolidSystem>;
     using SST = GetPropType<TypeTag, Properties::SolidState>;
     using MT = GetPropType<TypeTag, Properties::ModelTraits>;
-    using DTT = GetPropType<TypeTag, Properties::DispersionTensorType>;
     using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
     using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
 
-    using Traits = TracerVolumeVariablesTraits<PV, FSY, SSY, SST, MT, DTT, DT, EDM>;
+    using Traits = TracerVolumeVariablesTraits<PV, FSY, SSY, SST, MT, DT, EDM>;
 public:
     using type = TracerVolumeVariables<Traits>;
 };
diff --git a/dumux/porousmediumflow/tracer/volumevariables.hh b/dumux/porousmediumflow/tracer/volumevariables.hh
index 11958636579aca4dafa32171016522c00cf45b80..d70333917b333796042a149099662c23739a78b9 100644
--- a/dumux/porousmediumflow/tracer/volumevariables.hh
+++ b/dumux/porousmediumflow/tracer/volumevariables.hh
@@ -63,8 +63,6 @@ class TracerVolumeVariables
     static constexpr int numFluidComps = ParentType::numFluidComponents();
 
 public:
-    //! Export the dispersion tensor law
-    using DispersionTensorType = typename Traits::DispersionTensorType;
     //! Export the fluid system type
     using FluidSystem = typename Traits::FluidSystem;
     //! Export the solid state type
diff --git a/test/porousmediumflow/1pnc/dispersion/properties.hh b/test/porousmediumflow/1pnc/dispersion/properties.hh
index 160a623eeb488fde71c8180fbb5d8cf81fd07ef9..bf8c75367ace85b6541cddc75bb64dbd5a77c52f 100644
--- a/test/porousmediumflow/1pnc/dispersion/properties.hh
+++ b/test/porousmediumflow/1pnc/dispersion/properties.hh
@@ -103,10 +103,7 @@ template<class TypeTag>
 struct EnableThermalDispersion<TypeTag, TTag::DispersionTest> { static constexpr bool value = NONISOTHERMAL; };
 
 template<class TypeTag>
-struct DispersionFluxType<TypeTag, TTag::DispersionTest> { using type = OnePDiffusiveDispersionFlux<TypeTag>; };
-
-template<class TypeTag>
-struct DispersionTensorType<TypeTag, TTag::DispersionTest>
+struct CompositionalDispersionModel<TypeTag, TTag::DispersionTest>
 {
 #if SCHEIDEGGER
     using type = ScheideggersDispersionTensor<TypeTag>;
diff --git a/test/porousmediumflow/1pnc/dispersion/spatialparams.hh b/test/porousmediumflow/1pnc/dispersion/spatialparams.hh
index 521a34272192e426a7d7f2935acd723c1624046d..e8459422cad507649f92bb4ad838653fa9cf6240 100644
--- a/test/porousmediumflow/1pnc/dispersion/spatialparams.hh
+++ b/test/porousmediumflow/1pnc/dispersion/spatialparams.hh
@@ -25,6 +25,7 @@
 #ifndef DUMUX_1PNC_TEST_SPATIAL_PARAMS_HH
 #define DUMUX_1PNC_TEST_SPATIAL_PARAMS_HH
 
+#include <dune/common/exceptions.hh>
 #include <dumux/porousmediumflow/properties.hh>
 #include <dumux/porousmediumflow/fvspatialparams1p.hh>
 
@@ -62,6 +63,30 @@ public:
         alphaL_ = getParam<Scalar>("Problem.AlphaL");
         alphaT_ = getParam<Scalar>("Problem.AlphaT");
         dispersionTensorCoefficients_ = getParam<std::vector<Scalar>>("Problem.DispersionTensor");
+
+        if (dispersionTensorCoefficients_.size() > 1)
+        {
+            if (dispersionTensorCoefficients_.size() != (dimWorld*dimWorld))
+                DUNE_THROW(Dune::InvalidStateException, "For anisotropic dispersion tensors, please list all entries (dim x dim).");
+
+            int k = 0;
+            for (int i = 0; i < dimWorld; i++)
+            {
+                for (int j = 0; j < dimWorld; j++)
+                {
+                    dispersionTensor_[i][j] = dispersionTensorCoefficients_[k];
+                    k++;
+                }
+            }
+        }
+        else
+        {
+            if (dispersionTensorCoefficients_.size() != 1)
+                DUNE_THROW(Dune::InvalidStateException, "For isotropic dispersion tensors, please one scalar value.");
+
+            for (int i = 0; i < dimWorld; i++)
+                dispersionTensor_[i][i] = dispersionTensorCoefficients_[0];
+        }
     }
 
     /*!
@@ -93,7 +118,7 @@ public:
      *
      * \param globalPos The global position
      */
-    std::array<Scalar, 2> dispersionAlphas(const GlobalPosition& globalPos) const
+    std::array<Scalar, 2> dispersionAlphas(const GlobalPosition& globalPos, int phaseIdx = 0, int compIdx = 0) const
     { return { alphaL_, alphaT_ }; }
 
     /*!
@@ -101,35 +126,8 @@ public:
      *
      * \param globalPos The global position
      */
-    DimWorldMatrix dispersionTensor(const GlobalPosition& globalPos) const
-    {
-        DimWorldMatrix dispersionTensor(0.0);
-        if (dispersionTensorCoefficients_.size() > 1)
-        {
-            if (dispersionTensorCoefficients_.size() != (dimWorld*dimWorld))
-                DUNE_THROW(Dune::InvalidStateException, "For anisotropic dispersion tensors, please list all entries (dim x dim).");
-
-            int k = 0;
-            for (int i = 0; i < dimWorld; i++)
-            {
-                for (int j = 0; j < dimWorld; j++)
-                {
-                    dispersionTensor[i][j] = dispersionTensorCoefficients_[k];
-                    k++;
-                }
-            }
-        }
-        else
-        {
-            if (dispersionTensorCoefficients_.size() != 1)
-                DUNE_THROW(Dune::InvalidStateException, "For isotropic dispersion tensors, please one scalar value.");
-
-            for (int i = 0; i < dimWorld; i++)
-                dispersionTensor[i][i] = dispersionTensorCoefficients_[0];
-        }
-
-        return dispersionTensor;
-    }
+    const DimWorldMatrix &dispersionTensor(const GlobalPosition& globalPos, int phaseIdx = 0, int compIdx = 0) const
+    { return dispersionTensor_; }
 
 private:
     Scalar permeability_;
@@ -137,6 +135,7 @@ private:
     Scalar alphaL_;
     Scalar alphaT_;
     std::vector<Scalar> dispersionTensorCoefficients_;
+    DimWorldMatrix dispersionTensor_;
 };
 
 } // end namespace Dumux
diff --git a/test/porousmediumflow/tracer/constvel/spatialparams.hh b/test/porousmediumflow/tracer/constvel/spatialparams.hh
index 6ce643ddb77b6f61d8f8999a89805b0e8f5b2ed7..b9a2d97781a23089a8886aa4d620937539de1221 100644
--- a/test/porousmediumflow/tracer/constvel/spatialparams.hh
+++ b/test/porousmediumflow/tracer/constvel/spatialparams.hh
@@ -108,13 +108,14 @@ public:
                * elemVolVars[fvGeometry.scv(scvf.insideScvIdx())].extrusionFactor();
     }
 
-
     /*!
      * \brief Defines the dispersion tensor \f$\mathrm{[-]}\f$.
      *
      * \param globalPos The global position
      */
-    std::array<Scalar, 2> dispersionAlphas(const GlobalPosition& globalPos) const
+    std::array<Scalar, 2> dispersionAlphas(const GlobalPosition& globalPos,
+                                           const int phaseIdx = 0,
+                                           const int compIdx = 0) const
     { return { alphaL_, alphaT_ }; }
 
 private: