diff --git a/etc/caseDicts/postProcessing/graphs/graphLayerAverage b/etc/caseDicts/postProcessing/graphs/graphLayerAverage
index 23223ac58a..429089f587 100644
--- a/etc/caseDicts/postProcessing/graphs/graphLayerAverage
+++ b/etc/caseDicts/postProcessing/graphs/graphLayerAverage
@@ -10,18 +10,22 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-patches         (<patchNames>);  // Patches from which layers extrude
+patches         (<patchNames>);       // Patches from which layers extrude
 
-zones           ();              // Zones from which layers extrude
+zones           ();                   // Zones from which layers extrude
 
-axis            distance;        // The independent variable of the graph. Can
-                                 // be "x", "y", "z", "xyz" (all coordinates
-                                 // written out), or "distance" (from the start
-                                 // point).
+axis            distance;             // The independent variable of the graph.
+                                      // Can be "x", "y", "z", "xyz" (all
+                                      // coordinates written out), or
+                                      // "distance" (from the start point).
 
-symmetric       false;           // Are the layers symmetric about the centre?
+symmetric       false;                // Are the layers symmetric about the
+                                      // centre?
 
-fields          (<fieldsNames>); // Fields to plot
+fields          (<fieldsNames>);      // Fields to average
+
+//weightField   <weightFieldName>;    // Field or fields with which to weight
+//weightFields  (<weightFieldNames>); // the average
 
 #includeEtc "caseDicts/postProcessing/graphs/graphLayerAverage.cfg"
 
diff --git a/src/functionObjects/field/layerAverage/layerAverage.C b/src/functionObjects/field/layerAverage/layerAverage.C
index 7a33dee0bd..e42edb8a09 100644
--- a/src/functionObjects/field/layerAverage/layerAverage.C
+++ b/src/functionObjects/field/layerAverage/layerAverage.C
@@ -126,10 +126,10 @@ void Foam::functionObjects::layerAverage::calcLayers()
     }
 
     // Sum number of entries per layer
-    layerCount_ = sum(scalarField(mesh_.nCells(), 1));
+    layerVolume_ = sum<scalar>(mesh_.V());
 
     // Average the cell centres
-    layerCentre_ = sum(mesh_.cellCentres())/layerCount_;
+    layerCentre_ = sum<vector>(mesh_.V()*mesh_.C())/layerVolume_;
 
     // If symmetric, keep only half of the coordinates
     if (symmetric_)
@@ -139,6 +139,34 @@ void Foam::functionObjects::layerAverage::calcLayers()
 }
 
 
+Foam::tmp<Foam::VolInternalField<Foam::scalar>>
+Foam::functionObjects::layerAverage::weight() const
+{
+    if (weightFields_.empty())
+    {
+        return tmp<VolInternalField<scalar>>();
+    }
+
+    tmp<VolInternalField<scalar>> tresult =
+        VolInternalField<scalar>::New
+        (
+            "weight",
+            mesh_,
+            dimensionedScalar(dimless, scalar(1))
+        );
+
+    forAll(weightFields_, i)
+    {
+        const VolInternalField<scalar>& weightField =
+            mesh_.lookupObject<VolInternalField<scalar>>(weightFields_[i]);
+
+        tresult.ref() *= weightField;
+    }
+
+    return tresult;
+}
+
+
 template<>
 Foam::vector
 Foam::functionObjects::layerAverage::symmetricCoeff<Foam::vector>() const
@@ -245,6 +273,13 @@ bool Foam::functionObjects::layerAverage::read(const dictionary& dict)
 
     fields_ = dict.lookup<wordList>("fields");
 
+    weightFields_ =
+        dict.found("weightFields")
+      ? dict.lookup<wordList>("weightFields")
+      : dict.found("weightField")
+      ? wordList(1, dict.lookup<word>("weightField"))
+      : wordList();
+
     formatter_ = setWriter::New(dict.lookup("setFormat"), dict);
 
     calcLayers();
@@ -255,7 +290,9 @@ bool Foam::functionObjects::layerAverage::read(const dictionary& dict)
 
 Foam::wordList Foam::functionObjects::layerAverage::fields() const
 {
-    return fields_;
+    wordList result(fields_);
+    result.append(weightFields_);
+    return result;
 }
 
 
@@ -267,6 +304,13 @@ bool Foam::functionObjects::layerAverage::execute()
 
 bool Foam::functionObjects::layerAverage::write()
 {
+    // Get the weights
+    tmp<VolInternalField<scalar>> weight(this->weight());
+    tmp<Field<scalar>> layerWeight
+    (
+        weight.valid() ? sum<scalar>(mesh_.V()*weight) : tmp<Field<scalar>>()
+    );
+
     // Create list of available fields
     wordList fieldNames;
     forAll(fields_, fieldi)
@@ -275,7 +319,7 @@ bool Foam::functionObjects::layerAverage::write()
         (
             false
             #define FoundTypeField(Type, nullArg) \
-              || foundObject<VolField<Type>>(fields_[fieldi])
+              || foundObject<VolInternalField<Type>>(fields_[fieldi])
             FOR_ALL_FIELD_TYPES(FoundTypeField)
             #undef FoundTypeField
         )
@@ -295,16 +339,18 @@ bool Foam::functionObjects::layerAverage::write()
     #undef DeclareTypeValueSets
     forAll(fieldNames, fieldi)
     {
+        const word& fieldName = fieldNames[fieldi];
+
         #define CollapseTypeFields(Type, nullArg)                           \
-            if (mesh_.foundObject<VolField<Type>>(fieldNames[fieldi]))      \
+            if (mesh_.foundObject<VolInternalField<Type>>(fieldName))       \
             {                                                               \
-                const VolField<Type>& field =                               \
-                    mesh_.lookupObject<VolField<Type>>(fieldNames[fieldi]); \
+                const VolInternalField<Type>& field =                       \
+                    mesh_.lookupObject<VolInternalField<Type>>(fieldName);  \
                                                                             \
                 Type##ValueSets.set                                         \
                 (                                                           \
                     fieldi,                                                 \
-                    average(field.primitiveField())                         \
+                    average<Type>(weight, layerWeight, field)               \
                 );                                                          \
             }
         FOR_ALL_FIELD_TYPES(CollapseTypeFields);
diff --git a/src/functionObjects/field/layerAverage/layerAverage.H b/src/functionObjects/field/layerAverage/layerAverage.H
index 3634a961a7..17e69ea97c 100644
--- a/src/functionObjects/field/layerAverage/layerAverage.H
+++ b/src/functionObjects/field/layerAverage/layerAverage.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,7 +59,9 @@ Usage
         axis          | Component of the position to plot against | yes |
         symmetric     | Is the geometry symmetric around the centre layer? \
                                                    | no         | false
-        field         | Fields to average and plot | yes        |
+        fields        | Fields to average and plot | yes        |
+        weightField   | Field with which to weight the average | no | none
+        weightFields  | Fields with which to weight the average | no | ()
     \endtable
 
 SourceFiles
@@ -74,6 +76,7 @@ SourceFiles
 #include "fvMeshFunctionObject.H"
 #include "setWriter.H"
 #include "boolList.H"
+#include "volFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -113,15 +116,18 @@ class layerAverage
         //- Per cell the global layer
         labelList cellLayer_;
 
-        //- Per global layer the number of cells
-        scalarField layerCount_;
+        //- Per global layer the volume
+        scalarField layerVolume_;
 
         //- The average centre of each layer
         pointField layerCentre_;
 
-        //- Fields to sample
+        //- Fields to average
         wordList fields_;
 
+        //- Fields with which to weight the averages
+        wordList weightFields_;
+
         //- Set formatter
         autoPtr<setWriter> formatter_;
 
@@ -131,17 +137,26 @@ class layerAverage
         //- Create the layer information, the sort map, and the scalar axis
         void calcLayers();
 
+        //- Calculate and return the weight field, or a null pointer if there
+        //  are no weight fields
+        tmp<VolInternalField<scalar>> weight() const;
+
         //- Return the coefficient to multiply onto symmetric values
         template<class T>
         T symmetricCoeff() const;
 
         //- Sum field per layer
         template<class T>
-        Field<T> sum(const Field<T>& cellField) const;
+        tmp<Field<T>> sum(const VolInternalField<T>& cellField) const;
 
         //- Average a field per layer
         template<class T>
-        Field<T> average(const Field<T>& cellField) const;
+        tmp<Field<T>> average
+        (
+            const tmp<VolInternalField<scalar>>& cellWeight,
+            const tmp<Field<scalar>>& layerWeight,
+            const VolInternalField<T>& cellField
+        ) const;
 
 
 public:
diff --git a/src/functionObjects/field/layerAverage/layerAverageTemplates.C b/src/functionObjects/field/layerAverage/layerAverageTemplates.C
index 27154e3b62..18dfe2868a 100644
--- a/src/functionObjects/field/layerAverage/layerAverageTemplates.C
+++ b/src/functionObjects/field/layerAverage/layerAverageTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "layerAverage.H"
+#include "fvMesh.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -35,12 +36,13 @@ T Foam::functionObjects::layerAverage::symmetricCoeff() const
 
 
 template<class T>
-Foam::Field<T> Foam::functionObjects::layerAverage::sum
+Foam::tmp<Foam::Field<T>> Foam::functionObjects::layerAverage::sum
 (
-    const Field<T>& cellField
+    const VolInternalField<T>& cellField
 ) const
 {
-    Field<T> layerField(nLayers_, Zero);
+    tmp<Field<T>> tlayerField(new Field<T>(nLayers_, Zero));
+    Field<T>& layerField = tlayerField.ref();
 
     forAll(cellLayer_, celli)
     {
@@ -53,22 +55,30 @@ Foam::Field<T> Foam::functionObjects::layerAverage::sum
     Pstream::listCombineGather(layerField, plusEqOp<T>());
     Pstream::listCombineScatter(layerField);
 
-    return layerField;
+    return tlayerField;
 }
 
 
 template<class T>
-Foam::Field<T> Foam::functionObjects::layerAverage::average
+Foam::tmp<Foam::Field<T>> Foam::functionObjects::layerAverage::average
 (
-    const Field<T>& cellField
+    const tmp<VolInternalField<scalar>>& cellWeight,
+    const tmp<Field<scalar>>& layerWeight,
+    const VolInternalField<T>& cellField
 ) const
 {
-    // Sum and average
-    Field<T> layerField(sum(cellField)/layerCount_);
+    tmp<Field<T>> tlayerField
+    (
+        cellWeight.valid()
+      ? sum<T>(mesh_.V()*cellWeight*cellField)/layerWeight
+      : sum<T>(mesh_.V()*cellField)/layerVolume_
+    );
 
     // Handle symmetry
     if (symmetric_)
     {
+        Field<T>& layerField = tlayerField.ref();
+
         const T coeff = symmetricCoeff<T>();
 
         for (label i=0; i<nLayers_/2; i++)
@@ -82,7 +92,7 @@ Foam::Field<T> Foam::functionObjects::layerAverage::average
         layerField.setSize(nLayers_/2);
     }
 
-    return layerField;
+    return tlayerField;
 }