diff --git a/src/OpenFOAM/primitives/functions/Function1/Function1/Function1.H b/src/OpenFOAM/primitives/functions/Function1/Function1/Function1.H
index 121165fa05..33b10d3976 100644
--- a/src/OpenFOAM/primitives/functions/Function1/Function1/Function1.H
+++ b/src/OpenFOAM/primitives/functions/Function1/Function1/Function1.H
@@ -243,6 +243,9 @@ public:
         //- Is value constant (i.e. independent of x)
         virtual bool constant() const { return false; }
 
+        //- Can function be evaluated?
+        virtual bool good() const { return true; }
+
 
     // Evaluation
 
diff --git a/src/OpenFOAM/primitives/functions/Function1/None/NoneFunction1.H b/src/OpenFOAM/primitives/functions/Function1/None/NoneFunction1.H
index 64ffd6e428..2faa535cfb 100644
--- a/src/OpenFOAM/primitives/functions/Function1/None/NoneFunction1.H
+++ b/src/OpenFOAM/primitives/functions/Function1/None/NoneFunction1.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,8 +34,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Function1Types_None_H
-#define Function1Types_None_H
+#ifndef Foam_Function1Types_None_H
+#define Foam_Function1Types_None_H
 
 #include "Function1.H"
 
@@ -102,6 +102,9 @@ public:
         //- Value is independent of x
         virtual inline bool constant() const { return true; }
 
+        //- Function cannot be evaluated
+        virtual inline bool good() const { return false; }
+
         //- Placeholder: generates an error if called
         virtual Type value(const scalar x) const;
 
diff --git a/src/faOptions/faOption/faOption.C b/src/faOptions/faOption/faOption.C
index 89632a88c0..e95fd30c44 100644
--- a/src/faOptions/faOption/faOption.C
+++ b/src/faOptions/faOption/faOption.C
@@ -44,7 +44,7 @@ namespace Foam
 
 void Foam::fa::option::resetApplied()
 {
-    applied_.resize(fieldNames_.size());
+    applied_.resize_nocopy(fieldNames_.size());
     applied_ = false;
 }
 
@@ -91,7 +91,7 @@ Foam::autoPtr<Foam::fa::option> Foam::fa::option::New
     Info<< indent
         << "Selecting finite area options type " << modelType << endl;
 
-    const_cast<Time&>(patch.boundaryMesh().mesh().time()).libs().open
+    patch.boundaryMesh().mesh().time().libs().open
     (
         coeffs,
         "libs",
@@ -102,15 +102,16 @@ Foam::autoPtr<Foam::fa::option> Foam::fa::option::New
 
     if (!ctorPtr)
     {
-        FatalErrorInFunction
-            << "Unknown faOption model type "
-            << modelType << nl << nl
-            << "Valid faOption types are:" << nl
-            << dictionaryConstructorTablePtr_->sortedToc()
-            << exit(FatalError);
+        FatalIOErrorInLookup
+        (
+            coeffs,
+            "faOption",
+            modelType,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
     }
 
-    return autoPtr<option>(ctorPtr(name, modelType, coeffs, patch));
+    return autoPtr<fa::option>(ctorPtr(name, modelType, coeffs, patch));
 }
 
 
diff --git a/src/faOptions/faOption/faOption.H b/src/faOptions/faOption/faOption.H
index e7bddf348f..8e55e5c29a 100644
--- a/src/faOptions/faOption/faOption.H
+++ b/src/faOptions/faOption/faOption.H
@@ -72,8 +72,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef faOption_H
-#define faOption_H
+#ifndef Foam_faOption_H
+#define Foam_faOption_H
 
 #include "faMatricesFwd.H"
 #include "areaFieldsFwd.H"
@@ -263,7 +263,7 @@ public:
             //- Return dictionary
             inline const dictionary& coeffs() const noexcept;
 
-            //- Return const access to the source active flag
+            //- True if source is active
             inline bool active() const noexcept;
 
             //- Set the applied flag to true for field index fieldi
diff --git a/src/finiteVolume/cfdTools/general/fvOptions/fvOption.C b/src/finiteVolume/cfdTools/general/fvOptions/fvOption.C
index d58183d00b..c6813778c5 100644
--- a/src/finiteVolume/cfdTools/general/fvOptions/fvOption.C
+++ b/src/finiteVolume/cfdTools/general/fvOptions/fvOption.C
@@ -47,7 +47,7 @@ namespace Foam
 
 void Foam::fv::option::resetApplied()
 {
-    applied_.resize(fieldNames_.size());
+    applied_.resize_nocopy(fieldNames_.size());
     applied_ = false;
 }
 
diff --git a/src/finiteVolume/cfdTools/general/fvOptions/fvOption.H b/src/finiteVolume/cfdTools/general/fvOptions/fvOption.H
index f7a7408b84..c8bdef7909 100644
--- a/src/finiteVolume/cfdTools/general/fvOptions/fvOption.H
+++ b/src/finiteVolume/cfdTools/general/fvOptions/fvOption.H
@@ -68,14 +68,13 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef fvOption_H
-#define fvOption_H
+#ifndef Foam_fvOption_H
+#define Foam_fvOption_H
 
 #include "fvMatricesFwd.H"
 #include "primitiveFieldsFwd.H"
 #include "volFieldsFwd.H"
 #include "dictionary.H"
-#include "Switch.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -83,6 +82,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class fvMesh;
 
 namespace fv
diff --git a/src/functionObjects/solvers/scalarTransport/scalarTransport.H b/src/functionObjects/solvers/scalarTransport/scalarTransport.H
index e22cf8dc4b..25593a5ca2 100644
--- a/src/functionObjects/solvers/scalarTransport/scalarTransport.H
+++ b/src/functionObjects/solvers/scalarTransport/scalarTransport.H
@@ -100,14 +100,12 @@ Usage
                     type            scalarSemiImplicitSource;
                     enabled         true;
 
-                    scalarSemiImplicitSourceCoeffs
+                    selectionMode   all;
+                    volumeMode      specific;
+
+                    sources
                     {
-                        selectionMode   all;
-                        volumeMode      specific;
-                        injectionRateSuSp
-                        {
-                            s           (1 0);
-                        }
+                        s           (1 0);
                     }
                 }
             }
diff --git a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C
index 45698769a7..98b83b8800 100644
--- a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C
+++ b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.C
@@ -50,56 +50,159 @@ Foam::fv::SemiImplicitSource<Type>::volumeModeTypeNames_
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class Type>
-void Foam::fv::SemiImplicitSource<Type>::setFieldInjectionRates
+void Foam::fv::SemiImplicitSource<Type>::setFieldCoeffs
 (
-    const dictionary& dict
+    const dictionary& dict,
+    const word& keyExplicit,
+    const word& keyImplicit
 )
 {
     label count = dict.size();
 
     fieldNames_.resize_nocopy(count);
 
-    Su_.resize(count);
-    Sp_.resize(count);
+    Su_.clear();
+    Sp_.clear();
+    Su_.resize(2*count);
+    Sp_.resize(2*count);
+
+    driverSu_.clear();
+    driverSp_.clear();
+    driverSu_.resize(2*count);
+    driverSp_.resize(2*count);
+
+    valueExprSu_.clear();
+    valueExprSp_.clear();
+    valueExprSu_.resize(2*count);
+    valueExprSp_.resize(2*count);
 
     fv::option::resetApplied();
 
+    word modelType;
+    Tuple2<Type, scalar> sourceRates;
+
     count = 0;
+
     for (const entry& dEntry : dict)
     {
         const word& fieldName = dEntry.keyword();
+        bool ok = false;
 
         if (dEntry.isDict())
         {
-            const dictionary& subdict = dEntry.dict();
-            Su_.set(count, Function1<Type>::New("Su", subdict, &mesh_));
-            Sp_.set(count, Function1<scalar>::New("Sp", subdict, &mesh_));
+            const dictionary& subDict = dEntry.dict();
+
+            const entry* eptr;
+
+            if
+            (
+                (eptr = subDict.findEntry(keyExplicit, keyType::LITERAL))
+             != nullptr
+            )
+            {
+                ok = true;
+
+                if
+                (
+                    eptr->isDict()
+                 && eptr->dict().readEntry("type", modelType, keyType::LITERAL)
+                 && (modelType == "exprField")
+                )
+                {
+                    const dictionary& exprDict = eptr->dict();
+
+                    valueExprSu_.emplace_set(fieldName);
+                    valueExprSu_[fieldName].readEntry("expression", exprDict);
+
+                    driverSu_.set
+                    (
+                        fieldName,
+                        new expressions::volumeExprDriver(mesh_, exprDict)
+                    );
+                }
+                else
+                {
+                    Su_.set
+                    (
+                        fieldName,
+                        Function1<Type>::New(keyExplicit, subDict, &mesh_)
+                    );
+                }
+            }
+
+            if
+            (
+                (eptr = subDict.findEntry(keyImplicit, keyType::LITERAL))
+             != nullptr
+            )
+            {
+                ok = true;
+
+                if
+                (
+                    eptr->isDict()
+                 && eptr->dict().readEntry("type", modelType, keyType::LITERAL)
+                 && (modelType == "exprField")
+                )
+                {
+                    const dictionary& exprDict = eptr->dict();
+
+                    valueExprSp_.emplace_set(fieldName);
+                    valueExprSp_[fieldName].readEntry("expression", exprDict);
+
+                    driverSp_.set
+                    (
+                        fieldName,
+                        new expressions::volumeExprDriver(mesh_, exprDict)
+                    );
+                }
+                else
+                {
+                    Sp_.set
+                    (
+                        fieldName,
+                        Function1<scalar>::New(keyImplicit, subDict, &mesh_)
+                    );
+                }
+            }
         }
         else
         {
-            Tuple2<Type, scalar> injectionRate;
-            dEntry.readEntry(injectionRate);
+            // Non-dictionary form
+
+            dEntry.readEntry(sourceRates);
+
+            ok = true;
 
             Su_.set
             (
-                count,
+                fieldName,
                 new Function1Types::Constant<Type>
                 (
-                    "Su",
-                    injectionRate.first()
+                    keyExplicit,
+                    sourceRates.first()
                 )
             );
             Sp_.set
             (
-                count,
+                fieldName,
                 new Function1Types::Constant<scalar>
                 (
-                    "Sp",
-                    injectionRate.second()
+                    keyImplicit,
+                    sourceRates.second()
                 )
             );
         }
 
+        if (!ok)
+        {
+            FatalIOErrorInFunction(dict)
+                << "Require at least one of "
+                << keyExplicit << '/' << keyImplicit << " entries for "
+                << "field: " << fieldName << endl
+                << exit(FatalIOError);
+        }
+
         fieldNames_[count] = fieldName;
         ++count;
     }
@@ -154,6 +257,7 @@ void Foam::fv::SemiImplicitSource<Type>::addSup
 
     const GeometricField<Type, fvPatchField, volMesh>& psi = eqn.psi();
 
+    // Note: field name may deviate from psi name
     const word& fieldName = fieldNames_[fieldi];
 
     const scalar tmVal = mesh_.time().timeOutputValue();
@@ -161,67 +265,254 @@ void Foam::fv::SemiImplicitSource<Type>::addSup
     const dimensionSet SuDims(eqn.dimensions()/dimVolume);
     const dimensionSet SpDims(SuDims/psi.dimensions());
 
-
     // Explicit source
     {
-        const dimensioned<Type> SuValue
-        (
-            "Su",
-            SuDims,
-            Su_[fieldi].value(tmVal)/VDash_
-        );
+        const auto iter1 = valueExprSu_.cfind(fieldName);
+        const auto iter2 = Su_.cfind(fieldName);
 
-        if (mag(SuValue.value()) <= ROOTVSMALL)
+        tmp<DimensionedField<Type, volMesh>> tsu;
+
+        if (iter1.found())
         {
-            // No-op
-        }
-        else if (this->useSubMesh())
-        {
-            auto tsu = DimensionedField<Type, volMesh>::New
+            const auto& valueExpr = iter1.val();
+
+            typedef
+                GeometricField<Type, fvPatchField, volMesh>
+                ExprResultType;
+
+            if (debug)
+            {
+                Info<< "Explicit expression source:" << nl
+                    << ">>>>" << nl
+                    << valueExpr.c_str() << nl
+                    << "<<<<" << nl;
+            }
+
+            auto& driver = *(driverSu_[fieldName]);
+
+            driver.clearVariables();
+
+            if (notNull(rho))
+            {
+                driver.addContextObject("rho", &rho);
+            }
+
+            // Switch dimension checking off
+            const bool oldDimChecking = dimensionSet::checking(false);
+
+            driver.parse(valueExpr);
+
+            // Restore dimension checking
+            dimensionSet::checking(oldDimChecking);
+
+            const ExprResultType* ptr = driver.isResultType<ExprResultType>();
+
+            if (!ptr)
+            {
+                FatalErrorInFunction
+                    << "Expression for Su " << fieldName
+                    << " evaluated to <" << driver.resultType()
+                    << "> but expected <" << ExprResultType::typeName
+                    << ">" << endl
+                    << exit(FatalError);
+            }
+            else if (ptr->size() != mesh_.nCells())
+            {
+                FatalErrorInFunction
+                    << "Expression for Su " << fieldName
+                    << " evaluated to " << ptr->size()
+                    << " instead of " << mesh_.nCells() << " values" << endl
+                    << exit(FatalError);
+            }
+
+            if (notNull(rho))
+            {
+                driver.removeContextObject(&rho);
+            }
+
+            const Field<Type>& exprFld = ptr->primitiveField();
+
+            tsu = DimensionedField<Type, volMesh>::New
             (
                 name_ + fieldName + "Su",
                 mesh_,
                 dimensioned<Type>(SuDims, Zero)
             );
-            UIndirectList<Type>(tsu.ref(), cells_) = SuValue.value();
 
-            eqn += tsu;
+            if (this->useSubMesh())
+            {
+                for (const label celli : cells_)
+                {
+                    tsu.ref()[celli] = exprFld[celli]/VDash_;
+                }
+            }
+            else
+            {
+                tsu.ref().field() = exprFld;
+
+                if (!equal(VDash_, 1))
+                {
+                    tsu.ref().field() /= VDash_;
+                }
+            }
         }
-        else
+        else if (iter2.found() && iter2.val()->good())
         {
-            eqn += SuValue;
+            const dimensioned<Type> SuValue
+            (
+                "Su",
+                SuDims,
+                iter2.val()->value(tmVal)/VDash_
+            );
+
+            if (mag(SuValue.value()) <= ROOTVSMALL)
+            {
+                // No-op
+            }
+            else if (this->useSubMesh())
+            {
+                tsu = DimensionedField<Type, volMesh>::New
+                (
+                    name_ + fieldName + "Su",
+                    mesh_,
+                    dimensioned<Type>(SuDims, Zero)
+                );
+                UIndirectList<Type>(tsu.ref(), cells_) = SuValue.value();
+            }
+            else
+            {
+                eqn += SuValue;
+            }
+        }
+
+        if (tsu.valid())
+        {
+            eqn += tsu;
         }
     }
 
 
     // Implicit source
     {
-        const dimensioned<scalar> SpValue
-        (
-            "Sp",
-            SpDims,
-            Sp_[fieldi].value(tmVal)/VDash_
-        );
+        const auto iter1 = valueExprSp_.cfind(fieldName);
+        const auto iter2 = Sp_.cfind(fieldName);
 
-        if (mag(SpValue.value()) <= ROOTVSMALL)
+        tmp<DimensionedField<scalar, volMesh>> tsp;
+
+        if (iter1.found())
         {
-            // No-op
-        }
-        else if (this->useSubMesh())
-        {
-            auto tsp = DimensionedField<scalar, volMesh>::New
+            const auto& valueExpr = iter1.val();
+
+            typedef volScalarField ExprResultType;
+
+            if (debug)
+            {
+                Info<< "Implicit expression source:" << nl
+                    << ">>>>" << nl
+                    << valueExpr.c_str() << nl
+                    << "<<<<" << nl;
+            }
+
+            auto& driver = *(driverSp_[fieldName]);
+
+            driver.clearVariables();
+
+            if (notNull(rho))
+            {
+                driver.addContextObject("rho", &rho);
+            }
+
+            // Switch dimension checking off
+            const bool oldDimChecking = dimensionSet::checking(false);
+
+            driver.parse(valueExpr);
+
+            // Restore dimension checking
+            dimensionSet::checking(oldDimChecking);
+
+            const ExprResultType* ptr = driver.isResultType<ExprResultType>();
+
+            if (!ptr)
+            {
+                FatalErrorInFunction
+                    << "Expression for Sp " << fieldName
+                    << " evaluated to <" << driver.resultType()
+                    << "> but expected <" << ExprResultType::typeName
+                    << ">" << endl
+                    << exit(FatalError);
+            }
+            else if (ptr->size() != mesh_.nCells())
+            {
+                FatalErrorInFunction
+                    << "Expression for Sp " << fieldName
+                    << " evaluated to " << ptr->size()
+                    << " instead of " << mesh_.nCells() << " values" << endl
+                    << exit(FatalError);
+            }
+
+            if (notNull(rho))
+            {
+                driver.removeContextObject(&rho);
+            }
+
+            const Field<scalar>& exprFld = ptr->primitiveField();
+
+            tsp = DimensionedField<scalar, volMesh>::New
             (
                 name_ + fieldName + "Sp",
                 mesh_,
                 dimensioned<scalar>(SpDims, Zero)
             );
-            UIndirectList<scalar>(tsp.ref(), cells_) = SpValue.value();
 
-            eqn += fvm::SuSp(tsp, psi);
+            if (this->useSubMesh())
+            {
+                for (const label celli : cells_)
+                {
+                    tsp.ref()[celli] = exprFld[celli]/VDash_;
+                }
+            }
+            else
+            {
+                tsp.ref().field() = exprFld;
+
+                if (!equal(VDash_, 1))
+                {
+                    tsp.ref().field() /= VDash_;
+                }
+            }
         }
-        else
+        else if (iter2.found() && iter2.val()->good())
         {
-            eqn += fvm::SuSp(SpValue, psi);
+            const dimensioned<scalar> SpValue
+            (
+                "Sp",
+                SpDims,
+                iter2.val()->value(tmVal)/VDash_
+            );
+
+            if (mag(SpValue.value()) <= ROOTVSMALL)
+            {
+                // No-op
+            }
+            else if (this->useSubMesh())
+            {
+                tsp = DimensionedField<scalar, volMesh>::New
+                (
+                    name_ + fieldName + "Sp",
+                    mesh_,
+                    dimensioned<scalar>(SpDims, Zero)
+                );
+                UIndirectList<scalar>(tsp.ref(), cells_) = SpValue.value();
+            }
+            else
+            {
+                eqn += fvm::SuSp(SpValue, psi);
+            }
+        }
+
+        if (tsp.valid())
+        {
+            eqn += fvm::SuSp(tsp, psi);
         }
     }
 }
@@ -242,10 +533,26 @@ bool Foam::fv::SemiImplicitSource<Type>::read(const dictionary& dict)
             VDash_ = V_;
         }
 
+        // Compatibility (2112 and earlier)
+        const dictionary* injectDict =
+            coeffs_.findDict("injectionRateSuSp", keyType::LITERAL);
+
+        if (injectDict)
         {
-            setFieldInjectionRates
+            setFieldCoeffs
             (
-                coeffs_.subDict("injectionRateSuSp", keyType::LITERAL)
+                *injectDict,
+                "Su",  // Su = explicit
+                "Sp"   // Sp = implicit
+            );
+        }
+        else
+        {
+            setFieldCoeffs
+            (
+                coeffs_.subDict("sources", keyType::LITERAL),
+                "explicit",
+                "implicit"
             );
         }
 
diff --git a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H
index f51888d00c..fdcc83e75e 100644
--- a/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H
+++ b/src/fvOptions/sources/general/semiImplicitSource/SemiImplicitSource.H
@@ -33,7 +33,7 @@ Group
 Description
     Applies semi-implicit source within a specified region for \c Type,
     where \c \<Type\>=Scalar/Vector/SphericalTensor/SymmTensor/Tensor.
-    The injection rate coefficients are specified
+    The source rate coefficients are specified
     as pairs of Su-Sp coefficients, i.e.:
 
     \f[
@@ -57,27 +57,49 @@ Usage
 
         // Mandatory entries (runtime modifiable)
         volumeMode          <volumeModeType>;
-        injectionRateSuSp
+
+        // Specification of sources (OpenFOAM-2206 and newer)
+        sources
         {
+            // Specified as explicit(Su)/implicit(Sp) tuple:
             k           (30.7 0);
             epsilon     (1.5  0);
 
-            // The injectionRate can also be specified as a Function1
-            // by having dictionaries for the field entries instead
+            // Specified as Function1 or exprField
             k
             {
                 // Time-ramp from 0 to 30.7 at time 5
-                Su table
-                (
-                    (0 0.0)
-                    (5 30.7)
-                );
-                Sp      0.0;
+                explicit table ((0 0) (5 30.7));
+                implicit none;
             }
             epsilon
             {
-                Su      1.5;
-                Sp      0.0;
+                explicit
+                {
+                    type        exprField:
+                    expression  "(mag(pos()) < 0.1) ? 1.5 : 0";
+                }
+            }
+        }
+
+        // Traditional specification of sources (OpenFOAM-2112 and older)
+        injectionRateSuSp
+        {
+            // Specified as explicit(Su)/implicit(Sp) tuple:
+            k           (30.7 0);
+            epsilon     (1.5  0);
+
+            // Specified as Function1
+            k
+            {
+                // Time-ramp from 0 to 30.7 at time 5
+                Su  table ((0 0) (5 30.7));
+                Sp  0;
+            }
+            epsilon
+            {
+                Su  1.5;
+                Sp  0;
             }
         }
 
@@ -92,7 +114,8 @@ Usage
       type       | Type name: \<Type\>SemiImplicitSource <!--
                 -->                                     | word | yes  | -
       volumeMode | Volume mode type                     | word | yes  | -
-      injectionRateSuSp | Injection rate settings  | dictionary | yes | -
+      sources    | Explicit/implicit sources            | dict | cndtnl | -
+      injectionRateSuSp | Explicit/implicit sources     | dict | cndtnl | -
     \endtable
 
     The inherited entries are elaborated in:
@@ -105,6 +128,14 @@ Usage
       specific  | Values are given as \<quantity\>/m3
     \endverbatim
 
+Note
+    Missing explicit/implicit, Su/Sp entries are equivalent to constant values
+    of zero. However, at one entry must be supplied for the source terms.
+
+Note
+    To use the \c exprField (expression fields) handling, the \em sources
+    dictionary form must be used.
+
 See also
   - Foam::fvOption
 
@@ -119,6 +150,8 @@ SourceFiles
 #include "cellSetOption.H"
 #include "Enum.H"
 #include "Function1.H"
+#include "HashPtrTable.H"
+#include "volumeExprDriver.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -127,12 +160,6 @@ namespace Foam
 namespace fv
 {
 
-// Forward Declarations
-template<class Type> class SemiImplicitSource;
-
-template<class Type>
-Ostream& operator<<(Ostream&, const SemiImplicitSource<Type>&);
-
 /*---------------------------------------------------------------------------*\
                      Class SemiImplicitSource Declaration
 \*---------------------------------------------------------------------------*/
@@ -168,16 +195,34 @@ private:
         scalar VDash_;
 
         //- Explicit source contributions
-        PtrList<Function1<Type>> Su_;
+        HashPtrTable<Function1<Type>> Su_;
 
         //- Linearised implicit contributions
-        PtrList<Function1<scalar>> Sp_;
+        HashPtrTable<Function1<scalar>> Sp_;
+
+        //- Expression to evaluate for explicit source contribution
+        HashTable<expressions::exprString> valueExprSu_;
+
+        //- Expression to evaluate for linearised implicit contribution
+        HashTable<expressions::exprString> valueExprSp_;
+
+        //- Expression driver for explicit sources
+        HashPtrTable<expressions::volumeExprDriver> driverSu_;
+
+        //- Expression driver for implicit sources
+        HashPtrTable<expressions::volumeExprDriver> driverSp_;
 
 
     // Private Member Functions
 
-        //- Set the source coefficients from "injectionRateSuSp"
-        void setFieldInjectionRates(const dictionary& dict);
+        //- Set the source coefficients from "sources" (explicit/implicit)
+        //- or from "injectionRateSuSp" (Su/Sp)
+        void setFieldCoeffs
+        (
+            const dictionary& dict,
+            const word& keyExplicit,
+            const word& keyImplicit
+        );
 
 
 public:
diff --git a/src/meshTools/PatchFunction1/PatchFunction1/PatchFunction1.H b/src/meshTools/PatchFunction1/PatchFunction1/PatchFunction1.H
index 936b4512ce..2e3f6a6e98 100644
--- a/src/meshTools/PatchFunction1/PatchFunction1/PatchFunction1.H
+++ b/src/meshTools/PatchFunction1/PatchFunction1/PatchFunction1.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -218,6 +218,9 @@ public:
         //- Is value constant (i.e. independent of x)
         virtual bool constant() const { return false; }
 
+        //- Can function be evaluated?
+        virtual bool good() const { return true; }
+
         //- Is value uniform (i.e. independent of coordinate)
         virtual bool uniform() const = 0;
 
diff --git a/src/meshTools/PatchFunction1/UniformValueField/UniformValueField.H b/src/meshTools/PatchFunction1/UniformValueField/UniformValueField.H
index 2daafafd47..789aa3d70d 100644
--- a/src/meshTools/PatchFunction1/UniformValueField/UniformValueField.H
+++ b/src/meshTools/PatchFunction1/UniformValueField/UniformValueField.H
@@ -130,7 +130,13 @@ public:
         //- Is value constant (i.e. independent of x)
         virtual inline bool constant() const
         {
-            return uniformValuePtr_ && uniformValuePtr_->constant();
+            return (uniformValuePtr_ && uniformValuePtr_->constant());
+        }
+
+        //- Can function be evaluated?
+        virtual inline bool good() const
+        {
+            return (uniformValuePtr_ && uniformValuePtr_->good());
         }
 
         //- Is value uniform (i.e. independent of coordinate)
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/0.orig/tracer0 b/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/0.orig/tracer0
new file mode 100644
index 0000000000..f6224aa62d
--- /dev/null
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/0.orig/tracer0
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2206                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      tracer0;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    #includeEtc "caseDicts/setConstraintTypes"
+
+    rotor
+    {
+        type            zeroGradient;
+    }
+
+    stator
+    {
+        type            zeroGradient;
+    }
+
+    front
+    {
+        type            empty;
+    }
+
+    back
+    {
+        type            empty;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/scalarTransport b/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/scalarTransport
index 343a785f2b..c7dab79a6a 100644
--- a/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/scalarTransport
+++ b/tutorials/compressible/rhoPimpleFoam/RAS/mixerVessel2D/system/scalarTransport
@@ -33,9 +33,50 @@ tracer0
                 }
             }
 
-            injectionRateSuSp
+            sources
             {
-                tracer0     (0.01 0);
+                tracer0
+                {
+                    explicit constant 0.01;
+                }
+            }
+        }
+
+        source2
+        {
+            type            scalarSemiImplicitSource;
+
+            volumeMode      specific;
+
+            selectionMode   all;
+
+            sources
+            {
+                tracer0
+                {
+                    explicit
+                    {
+                        type       exprField;
+
+                        functions<scalar>
+                        {
+                            square
+                            {
+                                type square;
+                                scale 0.0025;
+                                level 0.0025;
+                                frequency 10;
+                            }
+                        }
+
+                        expression
+                        #{
+                            (hypot(pos().x() + 0.025, pos().y()) < 0.01)
+                          ? fn:square(time())
+                          : 0
+                        #};
+                    }
+                }
             }
         }
     }
diff --git a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/WatersKing/createFields.H b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/WatersKing/createFields.H
index dcc9f52656..f6d5cf49e0 100644
--- a/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/WatersKing/createFields.H
+++ b/tutorials/incompressible/pimpleFoam/laminar/planarPoiseuille/validation/WatersKing/createFields.H
@@ -49,10 +49,24 @@ IOdictionary fvOptions
         IOobject::NO_WRITE
     )
 );
-const dictionary& gradPDict =
-    fvOptions.subDict("momentumSource").subDict("injectionRateSuSp");
-const scalar K =
-    gradPDict.get<Tuple2<vector, scalar>>("U").first().x();
+
+scalar K(0);
+
+// Get x() component from U source
+{
+    const dictionary& momSource = fvOptions.subDict("momentumSource");
+
+    if (momSource.findDict("sources"))
+    {
+        K = momSource.subDict("sources")
+            .get<Tuple2<vector, scalar>>("U").first().x();
+    }
+    else
+    {
+        K = momSource.subDict("injectionRateSuSp")
+            .get<Tuple2<vector, scalar>>("U").first().x();
+    }
+}
 
 dictionary probes(IFstream(runTime.system()/"probes")());
 const point location = pointField(probes.lookup("probeLocations"))[0];