diff --git a/src/regionFaModels/thermalShell/thermalShell.C b/src/regionFaModels/thermalShell/thermalShell.C
index 797de7e2d5e90fc32f705d7e77f4887541a86be2..491b432be1fcb88821d230d91807e49f2a7919fe 100644
--- a/src/regionFaModels/thermalShell/thermalShell.C
+++ b/src/regionFaModels/thermalShell/thermalShell.C
@@ -43,16 +43,55 @@ defineTypeNameAndDebug(thermalShell, 0);
 
 addToRunTimeSelectionTable(thermalShellModel, thermalShell, dictionary);
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-bool thermalShell::read(const dictionary& dict)
+bool thermalShell::init(const dictionary& dict)
 {
+    if (thickness_ > 0)
+    {
+        h_ = dimensionedScalar("thickness", dimLength, thickness_);
+    }
+
     this->solution().readEntry("nNonOrthCorr", nNonOrthCorr_);
 
     return true;
 }
 
 
+tmp<areaScalarField> thermalShell::qr()
+{
+    IOobject io
+    (
+        "tqr",
+        primaryMesh().time().timeName(),
+        primaryMesh()
+    );
+
+    auto taqr =
+        tmp<areaScalarField>::New
+        (
+            io,
+            regionMesh(),
+            dimensionedScalar(dimPower/dimArea, Zero)
+        );
+
+    if (qrName_ != "none")
+    {
+        auto& aqr = taqr.ref();
+
+        const auto qr = primaryMesh().lookupObject<volScalarField>(qrName_);
+
+        const volScalarField::Boundary& vqr = qr.boundaryField();
+
+        aqr.primitiveFieldRef() = vsm().mapToSurface<scalar>(vqr);
+    }
+
+    return taqr;
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
 void thermalShell::solveEnergy()
 {
     if (debug)
@@ -68,7 +107,7 @@ void thermalShell::solveEnergy()
       - fam::laplacian(kappa()*h_, T_)
      ==
         qs_
-      //+ q(T_) handled by faOption contactHeatFlux
+      + qr()
       + faOptions()(h_, rhoCph, T_)
     );
 
@@ -118,18 +157,16 @@ thermalShell::thermalShell
             IOobject::AUTO_WRITE
         ),
         regionMesh()
-    )
+    ),
+    qrName_(dict.getOrDefault<word>("qr", "none")),
+    thickness_(dict.getOrDefault<scalar>("thickness", 0))
 {
-    init();
+    init(dict);
 }
 
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-void thermalShell::init()
-{}
-
-
 void thermalShell::preEvolveRegion()
 {}
 
diff --git a/src/regionFaModels/thermalShell/thermalShell.H b/src/regionFaModels/thermalShell/thermalShell.H
index 6160855c130e6b63aec44327d9f8f5ce421c0e5f..d675d68f3cf5c34be6847f1fdbbb6b0ac6bc0dae 100644
--- a/src/regionFaModels/thermalShell/thermalShell.H
+++ b/src/regionFaModels/thermalShell/thermalShell.H
@@ -98,7 +98,10 @@ class thermalShell
     // Private Member Functions
 
         //- Initialize thermalShell
-        void init();
+        bool init(const dictionary& dict);
+
+        //- Return radiative heat flux
+        tmp<areaScalarField> qr();
 
 
 protected:
@@ -125,12 +128,14 @@ protected:
             //- Thickness
             areaScalarField h_;
 
+            //- Name of the primary region radiative flux
+            const word qrName_;
 
-    // Protected Member Functions
+            //- Uniform thickness
+            scalar thickness_;
 
-        //- Read control parameters from dictionary
-        virtual bool read(const dictionary& dict);
 
+    // Protected Member Functions
 
         // Equations