diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
index 86499a7827..b7e52f3058 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.C
@@ -38,72 +38,9 @@ PopulationBalancePhaseSystem
     populationBalances_
     (
         this->lookup("populationBalances"),
-        diameterModels::populationBalanceModel::iNew(*this, dmdtfs_)
+        diameterModels::populationBalanceModel::iNew(*this)
     )
-{
-    forAll(populationBalances_, popBali)
-    {
-        const diameterModels::populationBalanceModel& popBal =
-            populationBalances_[popBali];
-
-        forAllConstIter
-        (
-            HashTable<const diameterModels::velocityGroup*>,
-            popBal.velocityGroupPtrs(),
-            iter1
-        )
-        {
-            const diameterModels::velocityGroup& velGrp1 = *iter1();
-
-            forAllConstIter
-            (
-                HashTable<const diameterModels::velocityGroup*>,
-                popBal.velocityGroupPtrs(),
-                iter2
-            )
-            {
-                const diameterModels::velocityGroup& velGrp2 = *iter2();
-
-                const phaseInterface interface
-                (
-                    velGrp1.phase(),
-                    velGrp2.phase()
-                );
-
-                if
-                (
-                    &velGrp1 != &velGrp2
-                    &&
-                    !dmdtfs_.found(interface)
-                )
-                {
-                    this->template validateMassTransfer
-                        <diameterModels::populationBalanceModel>(interface);
-
-                    dmdtfs_.insert
-                    (
-                        interface,
-                        new volScalarField
-                        (
-                            IOobject
-                            (
-                                IOobject::groupName
-                                (
-                                    "populationBalance:dmdtf",
-                                    interface.name()
-                                ),
-                                this->mesh().time().timeName(),
-                                this->mesh()
-                            ),
-                            this->mesh(),
-                            dimensionedScalar(dimDensity/dimTime, 0)
-                        )
-                    );
-                }
-            }
-        }
-    }
-}
+{}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -125,9 +62,12 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::dmdtf
 {
     tmp<volScalarField> tDmdtf = BasePhaseSystem::dmdtf(key);
 
-    if (dmdtfs_.found(key))
+    forAll(populationBalances_, popBali)
     {
-        tDmdtf.ref() += *dmdtfs_[key];
+        if (populationBalances_[popBali].dmdtfs().found(key))
+        {
+            tDmdtf.ref() += *populationBalances_[popBali].dmdtfs()[key];
+        }
     }
 
     return tDmdtf;
@@ -140,12 +80,20 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::dmdts() const
 {
     PtrList<volScalarField> dmdts(BasePhaseSystem::dmdts());
 
-    forAllConstIter(phaseSystem::dmdtfTable, dmdtfs_, dmdtfIter)
+    forAll(populationBalances_, popBali)
     {
-        const phaseInterface interface(*this, dmdtfIter.key());
+        forAllConstIter
+        (
+            phaseSystem::dmdtfTable,
+            populationBalances_[popBali].dmdtfs(),
+            dmdtfIter
+        )
+        {
+            const phaseInterface interface(*this, dmdtfIter.key());
 
-        addField(interface.phase1(), "dmdt", *dmdtfIter(), dmdts);
-        addField(interface.phase2(), "dmdt", - *dmdtfIter(), dmdts);
+            addField(interface.phase1(), "dmdt", *dmdtfIter(), dmdts);
+            addField(interface.phase2(), "dmdt", - *dmdtfIter(), dmdts);
+        }
     }
 
     return dmdts;
@@ -161,7 +109,10 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::momentumTransfer()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtUfs(dmdtfs_, eqns);
+    forAll(populationBalances_, popBali)
+    {
+        this->addDmdtUfs(populationBalances_[popBali].dmdtfs(), eqns);
+    }
 
     return eqnsPtr;
 }
@@ -176,7 +127,10 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::momentumTransferf()
 
     phaseSystem::momentumTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtUfs(dmdtfs_, eqns);
+    forAll(populationBalances_, popBali)
+    {
+        this->addDmdtUfs(populationBalances_[popBali].dmdtfs(), eqns);
+    }
 
     return eqnsPtr;
 }
@@ -191,7 +145,10 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::heatTransfer() const
 
     phaseSystem::heatTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtHefs(dmdtfs_, eqns);
+    forAll(populationBalances_, popBali)
+    {
+        this->addDmdtHefs(populationBalances_[popBali].dmdtfs(), eqns);
+    }
 
     return eqnsPtr;
 }
@@ -206,7 +163,10 @@ Foam::PopulationBalancePhaseSystem<BasePhaseSystem>::specieTransfer() const
 
     phaseSystem::specieTransferTable& eqns = eqnsPtr();
 
-    this->addDmdtYfs(dmdtfs_, eqns);
+    forAll(populationBalances_, popBali)
+    {
+        this->addDmdtYfs(populationBalances_[popBali].dmdtfs(), eqns);
+    }
 
     return eqnsPtr;
 }
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H
index 817128eb4e..e96150091a 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/PopulationBalancePhaseSystem/PopulationBalancePhaseSystem.H
@@ -62,9 +62,6 @@ class PopulationBalancePhaseSystem
         //- Population balances
         PtrList<diameterModels::populationBalanceModel> populationBalances_;
 
-        //- Mass transfer rates
-        phaseSystem::dmdtfTable dmdtfs_;
-
 
 public:
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
index 855d1c3441..d042116f1e 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H
@@ -181,10 +181,6 @@ protected:
             const phaseModelList& phaseModels
         ) const;
 
-        //- Check that mass transfer is supported across the given interface
-        template<class ModelType>
-        void validateMassTransfer(const phaseInterface& interface) const;
-
         //- Return the sum of the phase fractions of the moving phases
         tmp<volScalarField> sumAlphaMoving() const;
 
@@ -413,6 +409,10 @@ public:
             template<class ModelType>
             static const dictionary& modelSubDict(const dictionary& dict);
 
+            //- Check that mass transfer is supported across the given interface
+            template<class ModelType>
+            void validateMassTransfer(const phaseInterface& interface) const;
+
 
         // Sub-model lookup
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C
index 006b788117..9abbd002a1 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystemTemplates.C
@@ -28,23 +28,6 @@ License
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-template<class ModelType>
-void Foam::phaseSystem::validateMassTransfer
-(
-    const phaseInterface& interface
-) const
-{
-    if (interface.phase1().stationary() || interface.phase2().stationary())
-    {
-        FatalErrorInFunction
-            << "A " << ModelType::typeName << " was specified for pair "
-            << interface.name() << ", but one of these phases is stationary. "
-            << "Mass transfer is not supported on stationary phases"
-            << exit(FatalError);
-    }
-}
-
-
 namespace Foam
 {
 
@@ -455,6 +438,23 @@ const Foam::dictionary& Foam::phaseSystem::modelSubDict
 }
 
 
+template<class ModelType>
+void Foam::phaseSystem::validateMassTransfer
+(
+    const phaseInterface& interface
+) const
+{
+    if (interface.phase1().stationary() || interface.phase2().stationary())
+    {
+        FatalErrorInFunction
+            << "A " << ModelType::typeName << " was specified for pair "
+            << interface.name() << ", but one of these phases is stationary. "
+            << "Mass transfer is not supported on stationary phases"
+            << exit(FatalError);
+    }
+}
+
+
 template<class ModelType>
 bool Foam::phaseSystem::foundInterfacialModel
 (
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
index 53015d8555..6936ba60d7 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C
@@ -192,6 +192,63 @@ void Foam::diameterModels::populationBalanceModel::registerSizeGroups
 }
 
 
+void Foam::diameterModels::populationBalanceModel::initialiseDmdtfs()
+{
+    forAllConstIter
+    (
+        HashTable<const diameterModels::velocityGroup*>,
+        velocityGroupPtrs(),
+        iter1
+    )
+    {
+        const diameterModels::velocityGroup& velGrp1 = *iter1();
+
+        forAllConstIter
+        (
+            HashTable<const diameterModels::velocityGroup*>,
+            velocityGroupPtrs(),
+            iter2
+        )
+        {
+            const diameterModels::velocityGroup& velGrp2 = *iter2();
+
+            const phaseInterface interface(velGrp1.phase(), velGrp2.phase());
+
+            if
+            (
+                &velGrp1 != &velGrp2
+                &&
+                !dmdtfs_.found(interface)
+            )
+            {
+                fluid_.validateMassTransfer
+                    <diameterModels::populationBalanceModel>(interface);
+
+                dmdtfs_.insert
+                (
+                    interface,
+                    new volScalarField
+                    (
+                        IOobject
+                        (
+                            IOobject::groupName
+                            (
+                                "populationBalance:dmdtf",
+                                interface.name()
+                            ),
+                            mesh().time().timeName(),
+                            mesh()
+                        ),
+                        mesh(),
+                        dimensionedScalar(dimDensity/dimTime, 0)
+                    )
+                );
+            }
+        }
+    }
+}
+
+
 void Foam::diameterModels::populationBalanceModel::precompute()
 {
     forAll(coalescenceModels_, model)
@@ -260,22 +317,22 @@ void Foam::diameterModels::populationBalanceModel::birthByCoalescence
 
         const phaseInterface interfaceij(fi.phase(), fj.phase());
 
-        if (pDmdt_.found(interfaceij))
+        if (dmdtfs_.found(interfaceij))
         {
             const scalar dmdtSign =
                 interfaceij.index(fi.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interfaceij] += dmdtSign*fj.x()/v*Sui_*fj.phase().rho();
+            *dmdtfs_[interfaceij] += dmdtSign*fj.x()/v*Sui_*fj.phase().rho();
         }
 
         const phaseInterface interfaceik(fi.phase(), fk.phase());
 
-        if (pDmdt_.found(interfaceik))
+        if (dmdtfs_.found(interfaceik))
         {
             const scalar dmdtSign =
                 interfaceik.index(fi.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interfaceik] += dmdtSign*fk.x()/v*Sui_*fk.phase().rho();
+            *dmdtfs_[interfaceik] += dmdtSign*fk.x()/v*Sui_*fk.phase().rho();
         }
 
         sizeGroups_[i].shapeModelPtr()->addCoalescence(Sui_, fj, fk);
@@ -321,12 +378,12 @@ void Foam::diameterModels::populationBalanceModel::birthByBreakup
 
         const phaseInterface interface(fi.phase(), fk.phase());
 
-        if (pDmdt_.found(interface))
+        if (dmdtfs_.found(interface))
         {
             const scalar dmdtSign =
                 interface.index(fi.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interface] += dmdtSign*Sui_*fk.phase().rho();
+            *dmdtfs_[interface] += dmdtSign*Sui_*fk.phase().rho();
         }
 
         sizeGroups_[i].shapeModelPtr()->addBreakup(Sui_, fk);
@@ -398,12 +455,12 @@ void Foam::diameterModels::populationBalanceModel::birthByBinaryBreakup
 
     const phaseInterface interfaceij(fi.phase(), fj.phase());
 
-    if (pDmdt_.found(interfaceij))
+    if (dmdtfs_.found(interfaceij))
     {
         const scalar dmdtSign =
             interfaceij.index(fi.phase()) == 0 ? +1 : -1;
 
-        *pDmdt_[interfaceij] += dmdtSign*Sui_*fj.phase().rho();
+        *dmdtfs_[interfaceij] += dmdtSign*Sui_*fj.phase().rho();
     }
 
     dimensionedScalar Eta;
@@ -425,12 +482,12 @@ void Foam::diameterModels::populationBalanceModel::birthByBinaryBreakup
 
         const phaseInterface interfacekj(fk.phase(), fj.phase());
 
-        if (pDmdt_.found(interfacekj))
+        if (dmdtfs_.found(interfacekj))
         {
             const scalar dmdtSign =
                 interfacekj.index(fk.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interfacekj] += dmdtSign*Suk*fj.phase().rho();
+            *dmdtfs_[interfacekj] += dmdtSign*Suk*fj.phase().rho();
         }
 
         sizeGroups_[k].shapeModelPtr()->addBreakup(Suk, fj);
@@ -472,12 +529,12 @@ void Foam::diameterModels::populationBalanceModel::drift
 
         const phaseInterface interfacepe(fp.phase(), fe.phase());
 
-        if (pDmdt_.found(interfacepe))
+        if (dmdtfs_.found(interfacepe))
         {
             const scalar dmdtSign =
                 interfacepe.index(fp.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interfacepe] -= dmdtSign*Sue*fp.phase().rho();
+            *dmdtfs_[interfacepe] -= dmdtSign*Sue*fp.phase().rho();
         }
 
         sizeGroups_[i+1].shapeModelPtr()->addDrift(Sue, fp, model);
@@ -502,12 +559,12 @@ void Foam::diameterModels::populationBalanceModel::drift
 
         const phaseInterface interfacepw(fp.phase(), fw.phase());
 
-        if (pDmdt_.found(interfacepw))
+        if (dmdtfs_.found(interfacepw))
         {
             const scalar dmdtSign =
                 interfacepw.index(fp.phase()) == 0 ? +1 : -1;
 
-            *pDmdt_[interfacepw] -= dmdtSign*Suw*fp.phase().rho();
+            *dmdtfs_[interfacepw] -= dmdtSign*Suw*fp.phase().rho();
         }
 
         sizeGroups_[i-1].shapeModelPtr()->addDrift(Suw, fp, model);
@@ -547,7 +604,7 @@ void Foam::diameterModels::populationBalanceModel::sources()
         Sp_[i] = Zero;
     }
 
-    forAllIter(phaseSystem::dmdtfTable, pDmdt_, pDmdtIter)
+    forAllIter(phaseSystem::dmdtfTable, dmdtfs_, pDmdtIter)
     {
         *pDmdtIter() = Zero;
     }
@@ -720,8 +777,7 @@ bool Foam::diameterModels::populationBalanceModel::updateSources()
 Foam::diameterModels::populationBalanceModel::populationBalanceModel
 (
     const phaseSystem& fluid,
-    const word& name,
-    phaseSystem::dmdtfTable& pDmdt
+    const word& name
 )
 :
     regIOobject
@@ -734,7 +790,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
         )
     ),
     fluid_(fluid),
-    pDmdt_(pDmdt),
+    dmdtfs_(),
     mesh_(fluid_.mesh()),
     name_(name),
     dict_
@@ -812,6 +868,7 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
             << exit(FatalError);
     }
 
+    this->initialiseDmdtfs();
 
     if (coalescenceModels_.size() != 0)
     {
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
index 62ea1f1b83..b54b25d602 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H
@@ -240,7 +240,7 @@ private:
         const phaseSystem& fluid_;
 
         //- Interfacial mass transfer rates
-        phaseSystem::dmdtfTable& pDmdt_;
+        phaseSystem::dmdtfTable dmdtfs_;
 
         //- Reference to the mesh
         const fvMesh& mesh_;
@@ -336,6 +336,8 @@ private:
 
         void registerSizeGroups(sizeGroup& group);
 
+        void initialiseDmdtfs();
+
         void createPhasePairs();
 
         void precompute();
@@ -380,19 +382,10 @@ public:
     TypeName("populationBalanceModel");
 
 
-    // Constructor
+    // Constructors
 
-        populationBalanceModel
-        (
-            const phaseSystem& fluid,
-            const word& name,
-            HashPtrTable
-            <
-                volScalarField,
-                phaseInterfaceKey,
-                phaseInterfaceKey::hash
-            >& pDmdt
-        );
+        //- Construct for a fluid
+        populationBalanceModel(const phaseSystem& fluid, const word& name);
 
         //- Return clone
         autoPtr<populationBalanceModel> clone() const;
@@ -403,29 +396,22 @@ public:
         {
             const phaseSystem& fluid_;
 
-            phaseSystem::dmdtfTable& pDmdt_;
-
         public:
 
-            iNew
-            (
-                const phaseSystem& fluid,
-                phaseSystem::dmdtfTable& pDmdt
-            )
+            iNew(const phaseSystem& fluid)
             :
-                fluid_(fluid),
-                pDmdt_(pDmdt)
+                fluid_(fluid)
             {}
 
             autoPtr<populationBalanceModel> operator()(Istream& is) const
             {
-                word name(is);
+                const word name(is);
 
                 Info << "Setting up population balance: " << name << endl;
 
                 return autoPtr<populationBalanceModel>
                 (
-                    new populationBalanceModel(fluid_, name, pDmdt_)
+                    new populationBalanceModel(fluid_, name)
                 );
             }
         };
@@ -434,6 +420,7 @@ public:
     //- Destructor
     virtual ~populationBalanceModel();
 
+
     // Member Functions
 
         //- Dummy write for regIOobject
@@ -442,6 +429,9 @@ public:
         //- Return reference to the phaseSystem
         inline const phaseSystem& fluid() const;
 
+        //- Return reference to the interfacial mass transfer rates
+        inline const phaseSystem::dmdtfTable& dmdtfs() const;
+
         //- Return reference to the mesh
         inline const fvMesh& mesh() const;
 
diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
index 48b5c8be46..b87bdfa01c 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H
@@ -43,6 +43,13 @@ Foam::diameterModels::populationBalanceModel::fluid() const
 }
 
 
+inline const Foam::phaseSystem::dmdtfTable&
+Foam::diameterModels::populationBalanceModel::dmdtfs() const
+{
+    return dmdtfs_;
+}
+
+
 inline const Foam::fvMesh&
 Foam::diameterModels::populationBalanceModel::mesh() const
 {