diff --git a/src/regionModels/surfaceFilmModels/Make/files b/src/regionModels/surfaceFilmModels/Make/files
index 375ca36679..de4b1df758 100644
--- a/src/regionModels/surfaceFilmModels/Make/files
+++ b/src/regionModels/surfaceFilmModels/Make/files
@@ -47,6 +47,7 @@ $(THERMOMODELS)/phaseChangeModel/phaseChangeModel/phaseChangeModelNew.C
 $(THERMOMODELS)/phaseChangeModel/noPhaseChange/noPhaseChange.C
 $(THERMOMODELS)/phaseChangeModel/standardPhaseChange/standardPhaseChange.C
 $(THERMOMODELS)/phaseChangeModel/solidification/solidification.C
+$(THERMOMODELS)/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.C
 
 $(THERMOMODELS)/heatTransferModel/heatTransferModel/heatTransferModel.C
 $(THERMOMODELS)/heatTransferModel/heatTransferModel/heatTransferModelNew.C
@@ -67,6 +68,7 @@ $(THERMOMODELS)/filmViscosityModel/liquidViscosity/liquidViscosity.C
 $(THERMOMODELS)/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
 $(THERMOMODELS)/filmViscosityModel/ArrheniusViscosity/ArrheniusViscosity.C
 $(THERMOMODELS)/filmViscosityModel/function1Viscosity/function1Viscosity.C
+$(THERMOMODELS)/filmViscosityModel/waxSolventViscosity/waxSolventViscosity.C
 
 
 /* Boundary conditions */
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
index 548558cd26..b91b863829 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/thixotropicViscosity/thixotropicViscosity.C
@@ -154,7 +154,7 @@ void thixotropicViscosity::correct
         (
             max
             (
-                film.rhoSp(),
+               -film.rhoSp(),
                 dimensionedScalar("zero", film.rhoSp().dimensions(), 0)
             )/(deltaRho + deltaRho0),
             lambda_
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/waxSolventViscosity/waxSolventViscosity.C b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/waxSolventViscosity/waxSolventViscosity.C
new file mode 100644
index 0000000000..b3c4af2864
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/waxSolventViscosity/waxSolventViscosity.C
@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "waxSolventViscosity.H"
+#include "kinematicSingleLayer.H"
+#include "waxSolventEvaporation.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(waxSolventViscosity, 0);
+
+addToRunTimeSelectionTable
+(
+    filmViscosityModel,
+    waxSolventViscosity,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void waxSolventViscosity::correctMu()
+{
+    const kinematicSingleLayer& film = filmType<kinematicSingleLayer>();
+
+    const uniformDimensionedScalarField Wwax
+    (
+        film.regionMesh().lookupObject<uniformDimensionedScalarField>
+        (
+            waxSolventEvaporation::typeName + ":Wwax"
+        )
+    );
+
+    const uniformDimensionedScalarField Wsolvent
+    (
+        film.regionMesh().lookupObject<uniformDimensionedScalarField>
+        (
+            waxSolventEvaporation::typeName + ":Wsolvent"
+        )
+    );
+
+    const uniformDimensionedScalarField Ysolvent0
+    (
+        film.regionMesh().lookupObject<uniformDimensionedScalarField>
+        (
+            waxSolventEvaporation::typeName + ":Ysolvent0"
+        )
+    );
+
+    const volScalarField& Ysolvent
+    (
+        film.regionMesh().lookupObject<volScalarField>
+        (
+            waxSolventEvaporation::typeName + ":Ysolvent"
+        )
+    );
+
+    const volScalarField Xsolvent
+    (
+        Ysolvent*Wsolvent/((1 - Ysolvent)*Wwax + Ysolvent*Wsolvent)
+    );
+
+    const dimensionedScalar Xsolvent0
+    (
+        Ysolvent0*Wsolvent/((1 - Ysolvent0)*Wwax + Ysolvent0*Wsolvent)
+    );
+
+    mu_ = pow(muWax_/muSolvent_, (1 - Xsolvent)/(1 - Xsolvent0))*muSolvent_;
+    mu_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+waxSolventViscosity::waxSolventViscosity
+(
+    surfaceFilmRegionModel& film,
+    const dictionary& dict,
+    volScalarField& mu
+)
+:
+    filmViscosityModel(typeName, film, dict, mu),
+    muWax_
+    (
+        IOobject
+        (
+            typeName + ":muWax",
+            film.regionMesh().time().timeName(),
+            film.regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        film.regionMesh(),
+        dimensionedScalar("zero", dimDynamicViscosity, 0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    muWaxModel_
+    (
+        filmViscosityModel::New
+        (
+            film,
+            coeffDict_.subDict("muWax"),
+            muWax_
+        )
+    ),
+    muSolvent_
+    (
+        IOobject
+        (
+            typeName + ":muSolvent",
+            film.regionMesh().time().timeName(),
+            film.regionMesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        film.regionMesh(),
+        dimensionedScalar("zero", dimDynamicViscosity, 0),
+        zeroGradientFvPatchScalarField::typeName
+    ),
+    muSolventModel_
+    (
+        filmViscosityModel::New
+        (
+            film,
+            coeffDict_.subDict("muSolvent"),
+            muSolvent_
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+waxSolventViscosity::~waxSolventViscosity()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+void waxSolventViscosity::correct
+(
+    const volScalarField& p,
+    const volScalarField& T
+)
+{
+    muWaxModel_->correct(p, T);
+    muSolventModel_->correct(p, T);
+
+    correctMu();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/waxSolventViscosity/waxSolventViscosity.H b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/waxSolventViscosity/waxSolventViscosity.H
new file mode 100644
index 0000000000..8d82a294da
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/filmViscosityModel/waxSolventViscosity/waxSolventViscosity.H
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::regionModels::surfaceFilmModels::waxSolventViscosity
+
+Description
+    Wax solvent mixture viscosity model.
+
+SourceFiles
+    waxSolventViscosity.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef waxSolventViscosity_H
+#define waxSolventViscosity_H
+
+#include "filmViscosityModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class waxSolventViscosity Declaration
+\*---------------------------------------------------------------------------*/
+
+class waxSolventViscosity
+:
+    public filmViscosityModel
+{
+    // Private member functions
+
+        //- Correct the mixture viscosity
+        void correctMu();
+
+        //- Disallow default bitwise copy construct
+        waxSolventViscosity(const waxSolventViscosity&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const waxSolventViscosity&);
+
+
+protected:
+
+    // Protected data
+
+        //- Wax viscosity
+        volScalarField muWax_;
+
+        //- Wax viscosity model
+        autoPtr<filmViscosityModel> muWaxModel_;
+
+        //- Solvent viscosity
+        volScalarField muSolvent_;
+
+        //- Solvent viscosity model
+        autoPtr<filmViscosityModel> muSolventModel_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("waxSolvent");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        waxSolventViscosity
+        (
+            surfaceFilmRegionModel& film,
+            const dictionary& dict,
+            volScalarField& mu
+        );
+
+
+    //- Destructor
+    virtual ~waxSolventViscosity();
+
+
+    // Member Functions
+
+        //- Correct
+        virtual void correct
+        (
+            const volScalarField& p,
+            const volScalarField& T
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.C
new file mode 100644
index 0000000000..0228c4a54c
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.C
@@ -0,0 +1,406 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "waxSolventEvaporation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "thermoSingleLayer.H"
+#include "zeroField.H"
+
+#include "fvmDdt.H"
+#include "fvmDiv.H"
+#include "fvcDiv.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(waxSolventEvaporation, 0);
+
+addToRunTimeSelectionTable
+(
+    phaseChangeModel,
+    waxSolventEvaporation,
+    dictionary
+);
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+scalar waxSolventEvaporation::Sh
+(
+    const scalar Re,
+    const scalar Sc
+) const
+{
+    if (Re < 5.0E+05)
+    {
+        return 0.664*sqrt(Re)*cbrt(Sc);
+    }
+    else
+    {
+        return 0.037*pow(Re, 0.8)*cbrt(Sc);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+waxSolventEvaporation::waxSolventEvaporation
+(
+    surfaceFilmRegionModel& film,
+    const dictionary& dict
+)
+:
+    phaseChangeModel(typeName, film, dict),
+    Wwax_
+    (
+        IOobject
+        (
+            typeName + ":Wwax",
+            film.regionMesh().time().constant(),
+            film.regionMesh()
+        ),
+        readScalar(coeffDict_.lookup("Wwax"))
+    ),
+    Wsolvent_
+    (
+        IOobject
+        (
+            typeName + ":Wsolvent",
+            film.regionMesh().time().constant(),
+            film.regionMesh()
+        ),
+        readScalar(coeffDict_.lookup("Wsolvent"))
+    ),
+    Ysolvent0_
+    (
+        IOobject
+        (
+            typeName + ":Ysolvent0",
+            film.regionMesh().time().constant(),
+            film.regionMesh(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    Ysolvent_
+    (
+        IOobject
+        (
+            typeName + ":Ysolvent",
+            film.regionMesh().time().timeName(),
+            film.regionMesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        film.regionMesh()
+    ),
+    deltaMin_(readScalar(coeffDict_.lookup("deltaMin"))),
+    L_(readScalar(coeffDict_.lookup("L"))),
+    TbFactor_(coeffDict_.lookupOrDefault<scalar>("TbFactor", 1.1)),
+    YInfZero_(coeffDict_.lookupOrDefault<Switch>("YInfZero", false)),
+    activityCoeff_
+    (
+        Function1<scalar>::New("activityCoeff", coeffDict_)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+waxSolventEvaporation::~waxSolventEvaporation()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class YInfType>
+void waxSolventEvaporation::correctModel
+(
+    const scalar dt,
+    scalarField& availableMass,
+    scalarField& dMass,
+    scalarField& dEnergy,
+    YInfType YInf
+)
+{
+    const thermoSingleLayer& film = filmType<thermoSingleLayer>();
+
+    const volScalarField& delta = film.delta();
+    const volScalarField& deltaRho = film.deltaRho();
+    const surfaceScalarField& phi = film.phi();
+
+    // Set local thermo properties
+    const SLGThermo& thermo = film.thermo();
+    const filmThermoModel& filmThermo = film.filmThermo();
+    const label vapId = thermo.carrierId(filmThermo.name());
+
+    // Retrieve fields from film model
+    const scalarField& pInf = film.pPrimary();
+    const scalarField& T = film.T();
+    const scalarField& hs = film.hs();
+    const scalarField& rho = film.rho();
+    const scalarField& rhoInf = film.rhoPrimary();
+    const scalarField& muInf = film.muPrimary();
+    const scalarField& magSf = film.magSf();
+    const vectorField dU(film.UPrimary() - film.Us());
+    const scalarField limMass
+    (
+        max(scalar(0), availableMass - deltaMin_*rho*magSf)
+    );
+
+    // Molecular weight of vapour [kg/kmol]
+    const scalar Wvap = thermo.carrier().W(vapId);
+
+    const scalar Wwax = Wwax_.value();
+    const scalar Wsolvent = Wsolvent_.value();
+
+    volScalarField::Internal evapRateCoeff
+    (
+        IOobject
+        (
+            typeName + ":evapRateCoeff",
+            film.regionMesh().time().timeName(),
+            film.regionMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        film.regionMesh(),
+        dimensionedScalar("zero", dimDensity*dimVelocity, 0)
+    );
+
+    volScalarField::Internal evapRateInf
+    (
+        IOobject
+        (
+            typeName + ":evapRateInf",
+            film.regionMesh().time().timeName(),
+            film.regionMesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        film.regionMesh(),
+        dimensionedScalar("zero", dimDensity*dimVelocity, 0)
+    );
+
+    bool filmPresent = false;
+
+    forAll(dMass, celli)
+    {
+        if (delta[celli] > deltaMin_)
+        {
+            filmPresent = true;
+
+            const scalar Ysolvent = Ysolvent_[celli];
+
+            // Molefraction of solvent in liquid film
+            const scalar Xsolvent
+            (
+                Ysolvent*Wsolvent/((1 - Ysolvent)*Wwax + Ysolvent*Wsolvent)
+            );
+
+            // Primary region density [kg/m3]
+            const scalar rhoInfc = rhoInf[celli];
+
+            // Cell pressure [Pa]
+            const scalar pc = pInf[celli];
+
+            // Calculate the boiling temperature
+            const scalar Tb = filmThermo.Tb(pc);
+
+            // Local temperature - impose lower limit of 200 K for stability
+            const scalar Tloc = min(TbFactor_*Tb, max(200.0, T[celli]));
+
+            const scalar pPartialCoeff
+            (
+                filmThermo.pv(pc, Tloc)*activityCoeff_->value(Xsolvent)
+            );
+
+            scalar XsCoeff = pPartialCoeff/pc;
+
+            // Vapour phase mole fraction of solvent at interface
+            scalar Xs = XsCoeff*Xsolvent;
+
+            if (Xs > 1)
+            {
+                WarningInFunction
+                    << "Solvent vapour pressure > ambient pressure"
+                    << endl;
+
+                XsCoeff /= Xs;
+                Xs = 1;
+            }
+
+            // Vapour phase mass fraction of solvent at the interface
+            const scalar YsCoeff
+            (
+                XsCoeff/(XsCoeff*Xsolvent*Wsolvent + (1 - Xs)*Wvap)
+            );
+
+            // Primary region viscosity [Pa.s]
+            const scalar muInfc = muInf[celli];
+
+            // Reynolds number
+            const scalar Re = rhoInfc*mag(dU[celli])*L_/muInfc;
+
+            // Vapour diffusivity [m2/s]
+            const scalar Dab = filmThermo.D(pc, Tloc);
+
+            // Schmidt number
+            const scalar Sc = muInfc/(rhoInfc*(Dab + ROOTVSMALL));
+
+            // Sherwood number
+            const scalar Sh = this->Sh(Re, Sc);
+
+            // Mass transfer coefficient [m/s]
+            evapRateCoeff[celli] = rhoInfc*Sh*Dab/(L_ + ROOTVSMALL);
+
+            // Solvent mass transfer
+            const scalar dm
+            (
+                max
+                (
+                    dt*magSf[celli]
+                   *evapRateCoeff[celli]*(YsCoeff*Ysolvent - YInf[celli]),
+                    0
+                )
+            );
+
+            if (dm > limMass[celli])
+            {
+                evapRateCoeff[celli] *= limMass[celli]/dm;
+            }
+
+            evapRateInf[celli] = evapRateCoeff[celli]*YInf[celli];
+            evapRateCoeff[celli] *= YsCoeff;
+
+            // hVap[celli] = filmThermo.hl(pc, Tloc);
+        }
+    }
+
+    const dimensionedScalar deltaRho0Bydt
+    (
+        "deltaRho0",
+        deltaRho.dimensions()/dimTime,
+        ROOTVSMALL/dt
+    );
+
+    volScalarField::Internal impingementRate
+    (
+        max
+        (
+           -film.rhoSp()(),
+            dimensionedScalar("zero", film.rhoSp().dimensions(), 0)
+        )
+    );
+
+    if (filmPresent)
+    {
+        // Solve for the solvent mass fraction
+        fvScalarMatrix YsolventEqn
+        (
+            fvm::ddt(deltaRho, Ysolvent_)
+          + fvm::div(phi, Ysolvent_)
+         ==
+            deltaRho0Bydt*Ysolvent_()
+
+          + evapRateInf
+
+            // Include the effect of the impinging droplets
+            // added with Ysolvent = Ysolvent0
+          + impingementRate*Ysolvent0_
+
+          - fvm::Sp
+            (
+                deltaRho0Bydt
+              + evapRateCoeff
+              + film.rhoSp()()
+              + impingementRate,
+                Ysolvent_
+            )
+        );
+
+        YsolventEqn.relax();
+        YsolventEqn.solve();
+
+        Ysolvent_.min(1);
+        Ysolvent_.max(0);
+
+        scalarField dm
+        (
+            dt*magSf*rhoInf*(evapRateCoeff*Ysolvent_ + evapRateInf)
+        );
+
+        dMass += dm;
+
+        // Heat is assumed to be removed by heat-transfer to the wall
+        // so the energy remains unchanged by the phase-change.
+        dEnergy += dm*hs;
+
+        // Latent heat [J/kg]
+        // dEnergy += dm*(hs[celli] + hVap);
+    }
+}
+
+
+void waxSolventEvaporation::correctModel
+(
+    const scalar dt,
+    scalarField& availableMass,
+    scalarField& dMass,
+    scalarField& dEnergy
+)
+{
+    if (YInfZero_)
+    {
+        correctModel(dt, availableMass, dMass, dEnergy, zeroField());
+    }
+    else
+    {
+        const thermoSingleLayer& film = filmType<thermoSingleLayer>();
+        const label vapId = film.thermo().carrierId(film.filmThermo().name());
+        const scalarField& YInf = film.YPrimary()[vapId];
+
+        correctModel(dt, availableMass, dMass, dEnergy, YInf);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.H b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.H
new file mode 100644
index 0000000000..94d1a57e42
--- /dev/null
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.H
@@ -0,0 +1,160 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::waxSolventEvaporation
+
+Description
+    Wax solvent mixture evaporation model.
+
+SourceFiles
+    waxSolventEvaporation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef waxSolventEvaporation_H
+#define waxSolventEvaporation_H
+
+#include "phaseChangeModel.H"
+#include "uniformDimensionedFields.H"
+#include "Function1.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace regionModels
+{
+namespace surfaceFilmModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class waxSolventEvaporation Declaration
+\*---------------------------------------------------------------------------*/
+
+class waxSolventEvaporation
+:
+    public phaseChangeModel
+{
+    // Private member functions
+
+        //- Disallow default bitwise copy construct
+        waxSolventEvaporation(const waxSolventEvaporation&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const waxSolventEvaporation&);
+
+
+protected:
+
+    // Protected data
+
+        //- Molecular weight of wax [kg/kmol]
+        uniformDimensionedScalarField Wwax_;
+
+        //- Molecular weight of liquid [kg/kmol]
+        uniformDimensionedScalarField Wsolvent_;
+
+        //- Initial solvent mass-fraction
+        uniformDimensionedScalarField Ysolvent0_;
+
+        //- Solvent mass-fraction
+        volScalarField Ysolvent_;
+
+        //- Minimum film height for model to be active
+        const scalar deltaMin_;
+
+        //- Length scale [m]
+        const scalar L_;
+
+        //- Boiling temperature factor
+        //  Used to set max limit on temperature to Tb*TbFactor
+        const scalar TbFactor_;
+
+        //- Switch to treat YInf as zero
+        Switch YInfZero_;
+
+        //- Activity coefficient as a function of solvent mole fraction
+        autoPtr<Function1<scalar>> activityCoeff_;
+
+
+    // Protected member functions
+
+        //- Return Sherwood number as a function of Reynolds and Schmidt numbers
+        scalar Sh(const scalar Re, const scalar Sc) const;
+
+        template<class YInfType>
+        void correctModel
+        (
+            const scalar dt,
+            scalarField& availableMass,
+            scalarField& dMass,
+            scalarField& dEnergy,
+            YInfType YInf
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("waxSolventEvaporation");
+
+
+    // Constructors
+
+        //- Construct from surface film model
+        waxSolventEvaporation
+        (
+            surfaceFilmRegionModel& film,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~waxSolventEvaporation();
+
+
+    // Member Functions
+
+        //- Correct
+        virtual void correctModel
+        (
+            const scalar dt,
+            scalarField& availableMass,
+            scalarField& dMass,
+            scalarField& dEnergy
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFilmModels
+} // End namespace regionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
index 9f08a2a886..69ccdac22e 100644
--- a/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
+++ b/src/regionModels/surfaceFilmModels/thermoSingleLayer/thermoSingleLayer.C
@@ -628,6 +628,9 @@ void thermoSingleLayer::evolveRegion()
         InfoInFunction << endl;
     }
 
+    // Solve continuity for deltaRho_
+    solveContinuity();
+
     // Update sub-models to provide updated source contributions
     updateSubmodels();