diff --git a/src/atmosphericModels/Make/files b/src/atmosphericModels/Make/files
index 9e4f3018a5134d55d783acc2f4fcb9ddb2e89f26..62808d2debd74b05c7fe23be758d18ba691df0a0 100644
--- a/src/atmosphericModels/Make/files
+++ b/src/atmosphericModels/Make/files
@@ -33,6 +33,7 @@ fvOptions/atmPlantCanopyTurbSource/atmPlantCanopyTurbSource.C
fvOptions/atmPlantCanopyUSource/atmPlantCanopyUSource.C
fvOptions/atmPlantCanopyTSource/atmPlantCanopyTSource.C
fvOptions/atmNutSource/atmNutSource.C
+fvOptions/treeTurbulence/treeTurbulence.C
/* functionObjects */
diff --git a/src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.C b/src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.C
new file mode 100644
index 0000000000000000000000000000000000000000..6789d44f6e6f7194ef6a18ff2938b27dae9fd95f
--- /dev/null
+++ b/src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.C
@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "treeTurbulence.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+ defineTypeNameAndDebug(treeTurbulence, 0);
+ addToRunTimeSelectionTable(option, treeTurbulence, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
+
+Foam::volScalarField& Foam::fv::treeTurbulence::getOrReadField
+(
+ const word& fieldName
+) const
+{
+ auto* ptr = mesh_.getObjectPtr<volScalarField>(fieldName);
+
+ if (!ptr)
+ {
+ ptr = new volScalarField
+ (
+ IOobject
+ (
+ fieldName,
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh_
+ );
+ mesh_.objectRegistry::store(ptr);
+ }
+
+ return *ptr;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::fv::treeTurbulence::treeTurbulence
+(
+ const word& sourceName,
+ const word& modelType,
+ const dictionary& dict,
+ const fvMesh& mesh
+)
+:
+ fv::cellSetOption(sourceName, modelType, dict, mesh),
+ isEpsilon_(true),
+ betaP_(),
+ betaD_(),
+ Ceps1_(),
+ Ceps2_(),
+ betaStar_(),
+ CdName_(),
+ LADname_()
+{
+ read(dict);
+
+ const auto* turbPtr = mesh_.findObject<turbulenceModel>
+ (
+ turbulenceModel::propertiesName
+ );
+
+ if (!turbPtr)
+ {
+ FatalErrorInFunction
+ << "Unable to find a turbulence model."
+ << abort(FatalError);
+ }
+
+ fieldNames_.resize(2);
+
+ tmp<volScalarField> tepsilon = turbPtr->epsilon();
+ tmp<volScalarField> tomega = turbPtr->omega();
+
+ if (!tepsilon.isTmp())
+ {
+ fieldNames_[0] = tepsilon().name();
+ }
+ else if (!tomega.isTmp())
+ {
+ isEpsilon_ = false;
+ fieldNames_[0] = tomega().name();
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Unable to find epsilon or omega field." << nl
+ << abort(FatalError);
+ }
+
+ fieldNames_[1] = turbPtr->k()().name();
+
+ fv::option::resetApplied();
+
+ Log << " Applying treeTurbulence to: "
+ << fieldNames_[0] << " and " << fieldNames_[1]
+ << endl;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::fv::treeTurbulence::addSup
+(
+ fvScalarMatrix& eqn,
+ const label fieldi
+)
+{
+ if (fieldi == 1)
+ {
+ kSource(geometricOneField(), geometricOneField(), eqn, fieldi);
+ return;
+ }
+
+ if (isEpsilon_)
+ {
+ epsilonSource(geometricOneField(), geometricOneField(), eqn, fieldi);
+ }
+ else
+ {
+ omegaSource(geometricOneField(), geometricOneField(), eqn, fieldi);
+ }
+}
+
+
+void Foam::fv::treeTurbulence::addSup
+(
+ const volScalarField& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+)
+{
+ if (fieldi == 1)
+ {
+ kSource(geometricOneField(), rho, eqn, fieldi);
+ return;
+ }
+
+ if (isEpsilon_)
+ {
+ epsilonSource(geometricOneField(), rho, eqn, fieldi);
+ }
+ else
+ {
+ omegaSource(geometricOneField(), rho, eqn, fieldi);
+ }
+}
+
+
+void Foam::fv::treeTurbulence::addSup
+(
+ const volScalarField& alpha,
+ const volScalarField& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+)
+{
+ if (fieldi == 1)
+ {
+ kSource(alpha, rho, eqn, fieldi);
+ return;
+ }
+
+ if (isEpsilon_)
+ {
+ epsilonSource(alpha, rho, eqn, fieldi);
+ }
+ else
+ {
+ omegaSource(alpha, rho, eqn, fieldi);
+ }
+}
+
+
+bool Foam::fv::treeTurbulence::read(const dictionary& dict)
+{
+ if (!fv::cellSetOption::read(dict))
+ {
+ return false;
+ }
+
+ betaP_ = dict.getOrDefault<scalar>("betaP", 1.0);
+ betaD_ = dict.getOrDefault<scalar>("betaD", 4.0);
+ Ceps1_ = dict.getOrDefault<scalar>("Ceps1", 0.9);
+ Ceps2_ = dict.getOrDefault<scalar>("Ceps2", 0.9);
+ betaStar_ = dict.getOrDefault<scalar>("betaStar", 0.09);
+ CdName_ = dict.getOrDefault<word>("Cd", "Cd");
+ LADname_ = dict.getOrDefault<word>("LAD", "LAD");
+
+ (void) getOrReadField(CdName_);
+ (void) getOrReadField(LADname_);
+
+ return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.H b/src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.H
new file mode 100644
index 0000000000000000000000000000000000000000..478d2fc6d17e846904c97d53e73e617ec0324b3d
--- /dev/null
+++ b/src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulence.H
@@ -0,0 +1,262 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::fv::treeTurbulence
+
+Group
+ grpFvOptionsSources
+
+Description
+ Applies sources on \c k and either \c epsilon or \c omega
+ to incorporate effects of trees on atmospheric boundary layer.
+
+ Corrections applied to:
+ \verbatim
+ k | Turbulent kinetic energy [m^2/s^2]
+ \endverbatim
+
+ Corrections applied to either of the below, if exist:
+ \verbatim
+ epsilon | Turbulent kinetic energy dissipation rate [m^2/s^3]
+ omega | Specific dissipation rate [1/s]
+ \endverbatim
+
+ Required fields:
+ \verbatim
+ k | Turbulent kinetic energy [m^2/s^2]
+ Cd | Canopy drag coefficient [-]
+ LAD | Leaf area density [m^2/m^3]
+ epsilon | Turbulent kinetic energy dissipation rate [m^2/s^3]
+ omega | Specific dissipation rate [1/s]
+ \endverbatim
+
+ References:
+ \verbatim
+ Governing equations (tag:BSG):
+ Brozovsky, J., Simonsen, A., & Gaitani, N. (2021).
+ Validation of a CFD model for the evaluation of urban microclimate
+ at high latitudes: A case study in Trondheim, Norway.
+ Building and Environment, 205, 108175.
+ DOI:10.1016/j.buildenv.2021.108175
+ \endverbatim
+
+Usage
+ Example by using \c constant/fvOptions:
+ \verbatim
+ treeTurbulence1
+ {
+ // Mandatory entries
+ type treeTurbulence;
+
+ // Optional entries
+ Cd <word>;
+ LAD <word>;
+ betaP <scalar>;
+ betaD <scalar>;
+ Ceps1 <scalar>;
+ Ceps2 <scalar>;
+ betaStar <scalar>;
+
+ // Inherited entries
+ ...
+ }
+ \endverbatim
+
+ where the entries mean:
+ \table
+ Property | Description | Type | Reqd | Deflt
+ type | Type name: treeTurbulence | word | yes | -
+ Cd | Name of operand canopy drag coefficient field | word | no | Cd
+ LAD | Name of operand leaf area density field | word | no | LAD
+ betaP | Model coefficient | scalar | no | 1.0
+ betaD | Model coefficient | scalar | no | 4.0
+ Ceps1 | Model coefficient | scalar | no | 0.9
+ Ceps2 | Model coefficient | scalar | no | 0.9
+ betaStar | Model coefficient | scalar | no | 0.09
+ \endtable
+
+ The inherited entries are elaborated in:
+ - \link fvOption.H \endlink
+ - \link cellSetOption.H \endlink
+
+SourceFiles
+ treeTurbulence.C
+ treeTurbulenceTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fv_treeTurbulence_H
+#define fv_treeTurbulence_H
+
+#include "cellSetOption.H"
+#include "turbulentTransportModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+ Class treeTurbulence Declaration
+\*---------------------------------------------------------------------------*/
+
+class treeTurbulence
+:
+ public fv::cellSetOption
+{
+ // Private Data
+
+ //- Internal flag to determine the working field is epsilon or omega
+ bool isEpsilon_;
+
+ // Model coefficients
+
+ scalar betaP_;
+ scalar betaD_;
+ scalar Ceps1_;
+ scalar Ceps2_;
+ scalar betaStar_;
+
+ //- Name of operand canopy drag coefficient field
+ word CdName_;
+
+ //- Name of operand leaf area density field
+ word LADname_;
+
+
+ // Private Member Functions
+
+ //- Return requested field from the object registry
+ //- or read+register the field to the object registry
+ volScalarField& getOrReadField(const word& fieldName) const;
+
+ //- Apply sources to turbulent kinetic energy
+ template<class AlphaFieldType, class RhoFieldType>
+ void kSource
+ (
+ const AlphaFieldType& alpha,
+ const RhoFieldType& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+ ) const;
+
+ //- Apply sources to turbulent kinetic energy dissipation rate
+ template<class AlphaFieldType, class RhoFieldType>
+ void epsilonSource
+ (
+ const AlphaFieldType& alpha,
+ const RhoFieldType& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+ ) const;
+
+ //- Apply sources to specific dissipation rate
+ template<class AlphaFieldType, class RhoFieldType>
+ void omegaSource
+ (
+ const AlphaFieldType& alpha,
+ const RhoFieldType& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("treeTurbulence");
+
+
+ // Constructors
+
+ //- Construct from explicit source name and mesh
+ treeTurbulence
+ (
+ const word& sourceName,
+ const word& modelType,
+ const dictionary& dict,
+ const fvMesh& mesh
+ );
+
+ //- No copy construct
+ treeTurbulence(const treeTurbulence&) = delete;
+
+ //- No copy assignment
+ void operator=(const treeTurbulence&) = delete;
+
+
+ // Member Functions
+
+ //- Add explicit contribution to epsilon or omega equation
+ //- for incompressible flow computations
+ virtual void addSup
+ (
+ fvScalarMatrix& eqn,
+ const label fieldi
+ );
+
+ //- Add explicit contribution to epsilon or omega equation
+ //- for compressible flow computations
+ virtual void addSup
+ (
+ const volScalarField& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+ );
+
+ //- Add explicit contribution to epsilon or omega equation
+ //- for multiphase flow computations
+ virtual void addSup
+ (
+ const volScalarField& alpha,
+ const volScalarField& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+ );
+
+ //- Read source dictionary
+ virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+ #include "treeTurbulenceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulenceTemplates.C b/src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulenceTemplates.C
new file mode 100644
index 0000000000000000000000000000000000000000..58b6f6d477f44d9769a3e1a1832152e9fe2c3ef1
--- /dev/null
+++ b/src/atmosphericModels/fvOptions/treeTurbulence/treeTurbulenceTemplates.C
@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "treeTurbulence.H"
+#include "volFields.H"
+#include "fvmSup.H"
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::treeTurbulence::kSource
+(
+ const AlphaFieldType& alpha,
+ const RhoFieldType& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+) const
+{
+ const auto* turbPtr = mesh_.findObject<turbulenceModel>
+ (
+ turbulenceModel::propertiesName
+ );
+
+ const volScalarField& k = turbPtr->k();
+ const volVectorField& U = turbPtr->U();
+ const volScalarField magU(mag(U));
+ const volScalarField& Cd = getOrReadField(CdName_);
+ const volScalarField& LAD = getOrReadField(LADname_);
+
+ // (BSG:Eq. 8)
+ eqn +=
+ alpha*rho*LAD*Cd*betaP_*pow3(magU)
+ - fvm::Sp(alpha*rho*LAD*Cd*betaD_*magU, k);
+}
+
+
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::treeTurbulence::epsilonSource
+(
+ const AlphaFieldType& alpha,
+ const RhoFieldType& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+) const
+{
+ const auto* turbPtr = mesh_.findObject<turbulenceModel>
+ (
+ turbulenceModel::propertiesName
+ );
+
+ const volScalarField& epsilon = turbPtr->epsilon();
+ const volScalarField& k = turbPtr->k();
+ const volVectorField& U = turbPtr->U();
+ const volScalarField magU(mag(U));
+ const volScalarField& Cd = getOrReadField(CdName_);
+ const volScalarField& LAD = getOrReadField(LADname_);
+
+ // (BSG:Eq. 9)
+ eqn +=
+ fvm::Sp
+ (
+ alpha*rho*LAD*Cd*(Ceps1_*betaP_/k*pow3(magU) - Ceps2_*betaD_*magU),
+ epsilon
+ );
+}
+
+
+template<class AlphaFieldType, class RhoFieldType>
+void Foam::fv::treeTurbulence::omegaSource
+(
+ const AlphaFieldType& alpha,
+ const RhoFieldType& rho,
+ fvScalarMatrix& eqn,
+ const label fieldi
+) const
+{
+ const auto* turbPtr = mesh_.findObject<turbulenceModel>
+ (
+ turbulenceModel::propertiesName
+ );
+
+ const volScalarField& omega = turbPtr->omega();
+ const volScalarField& k = turbPtr->k();
+ const volVectorField& U = turbPtr->U();
+ const volScalarField magU(mag(U));
+ const volScalarField& Cd = getOrReadField(CdName_);
+ const volScalarField& LAD = getOrReadField(LADname_);
+
+ // (derived from BSG:Eq. 9 by assuming epsilon = betaStar*omega*k)
+ eqn +=
+ fvm::Sp
+ (
+ alpha*rho*LAD*Cd*betaStar_
+ *(Ceps1_*betaP_*pow3(magU) - Ceps2_*betaD_*k*magU),
+ omega
+ );
+}
+
+
+// ************************************************************************* //