diff --git a/CHANGELOG b/CHANGELOG
index c45f2392f29df6d7ed813dbeb16cc148430c9563..78f30cc31bb448e75b4901e61696a2150ea4f2eb 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -5,7 +5,9 @@ Differences Between DuMuX 2.8 and DuMuX 2.9
   -
 
 * IMPROVEMENTS and ENHANCEMENTS:
-  -
+  - The multidomain models have been restructured. This aims at a reduction in
+    duplicated code portions, more consistent procedure in the isothermal
+    and non.isothermal models, reduction of split functions.
 
 * IMMEDIATE INTERFACE CHANGES not allowing/requiring a deprecation period:
   - For the multidomain models, the notation of the boundary condition types
diff --git a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
index 78200c2cae3966517c550c46833df490f65d6274..0e90f01353fabf8aa17eada76c0d7644d011e670 100644
--- a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
+++ b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
@@ -220,7 +220,6 @@ public:
                                  cParams,
                                  couplingRes1, couplingRes2);
 
-        const GlobalPosition& globalPos1 = cParams.fvGeometry1.subContVol[vertInElem1].global;
         const GlobalPosition& globalPos2 = cParams.fvGeometry2.subContVol[vertInElem2].global;
 
         // ENERGY Balance
@@ -250,9 +249,6 @@ public:
                                         this->globalProblem().localResidual2().residual(vertInElem2)[energyEqIdx2]);
             }
         }
-
-        // TODO: unify the behavior for cParams.boundaryTypes2.isCouplingNeumann()
-        //       with the different one in the isothermal LOP
         if (cParams.boundaryTypes2.isCouplingNeumann(energyEqIdx2))
         {
             const GlobalPosition& bfNormal1 = boundaryVars1.face().normal;
@@ -293,7 +289,6 @@ public:
                                           * boundaryVars1.componentEnthalpy(phaseCompIdx1)
                                           * FluidSystem::molarMass(phaseCompIdx1);
 
-                // TODO: use mass transfer coefficient here?
                 couplingRes2.accumulate(lfsu2.child(energyEqIdx2), vertInElem2,
                                         -(convectiveFlux - sumDiffusiveEnergyFlux - conductiveFlux));
             }
diff --git a/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh b/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
index d2f8e4be1cf76cdad3fe5ca0d6cc3af681e3591a..12ceb922dbefeb1258bc72a48e41e60279bafd44 100644
--- a/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
+++ b/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
@@ -548,6 +548,9 @@ public:
             // only enter here, if a boundary layer model is used for the computation of the diffusive fluxes
             if (blModel_)
             {
+                Scalar advectiveFlux = normalMassFlux1
+                                       * cParams.elemVolVarsCur1[vertInElem1].massFraction(transportCompIdx1);
+
                 Scalar diffusiveFlux = bfNormal1.two_norm()
                                        * evalBoundaryLayerConcentrationGradient(cParams, vertInElem1)
                                        * (boundaryVars1.diffusionCoeff(transportCompIdx1)
@@ -555,12 +558,9 @@ public:
                                        * boundaryVars1.molarDensity()
                                        * FluidSystem::molarMass(transportCompIdx1);
 
-                Scalar advectiveFlux = normalMassFlux1
-                                       * cParams.elemVolVarsCur1[vertInElem1].massFraction(transportCompIdx1);
-
                 const Scalar massTransferCoeff = evalMassTransferCoefficient(cParams, vertInElem1, vertInElem2);
-                // TODO: unify this behavior with the one in the non-isothermal LOP
-                if (globalProblem_.sdProblem1().isCornerPoint(globalPos1) && massTransferModel_)
+
+                if (massTransferModel_ && globalProblem_.sdProblem1().isCornerPoint(globalPos1))
                 {
                     Scalar diffusiveFluxAtCorner = bfNormal1
                                                    * boundaryVars1.moleFractionGrad(transportCompIdx1)
@@ -587,6 +587,7 @@ public:
 
                 Scalar advectiveFlux = normalMassFlux1
                                        * cParams.elemVolVarsCur1[vertInElem1].massFraction(transportCompIdx1);
+
                 Scalar diffusiveFlux = bfNormal1
                                        * boundaryVars1.moleFractionGrad(transportCompIdx1)
                                        * (boundaryVars1.diffusionCoeff(transportCompIdx1)