diff --git a/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh
index 715b87cc97dcb67bd8f5fb77a86a7bf52716825e..3939f1f1cdca93fa4d535190350c566c611df088 100644
--- a/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh
+++ b/dumux/discretization/cellcentered/mpfa/localassemblerbase.hh
@@ -156,7 +156,7 @@ class InteractionVolumeAssemblerBase
             // gravitational acceleration on this face
             const auto& curLocalScvf = iv.localScvf(faceIdx);
             const auto& curGlobalScvf = fvGeometry().scvf(curLocalScvf.gridScvfIndex());
-            const auto& gravity = problem().gravityAtPos(curGlobalScvf.ipGlobal());
+            const auto& gravity = problem().spatialParams().gravity(curGlobalScvf.ipGlobal());
 
             // get permeability tensor in "positive" sub volume
             const auto& neighborScvIndices = curLocalScvf.neighboringLocalScvIndices();
diff --git a/dumux/flux/box/darcyslaw.hh b/dumux/flux/box/darcyslaw.hh
index 4abfcf0ecd799064c1f7cce605426ed5e06e9cc1..973fd3e8557da298c1663180ec0d7f882d06aaa1 100644
--- a/dumux/flux/box/darcyslaw.hh
+++ b/dumux/flux/box/darcyslaw.hh
@@ -117,7 +117,7 @@ public:
         }
 
         if (enableGravity)
-            gradP.axpy(-rho, problem.gravityAtPos(scvf.center()));
+            gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center()));
 
         // apply the permeability and return the flux
         return -1.0*vtmv(scvf.unitOuterNormal(), K, gradP)*scvf.area();
diff --git a/dumux/flux/cctpfa/darcyslaw.hh b/dumux/flux/cctpfa/darcyslaw.hh
index ab5843a73d9335768521d8ed2d78dda7e16dfc2c..41c1fc92fb2e3e818aa97e36fa45a7d13904ffe5 100644
--- a/dumux/flux/cctpfa/darcyslaw.hh
+++ b/dumux/flux/cctpfa/darcyslaw.hh
@@ -179,7 +179,7 @@ class CCTpfaDarcysLaw<ScalarType, FVGridGeometry, /*isNetwork*/ false>
             const auto pOutside = outsideVolVars.pressure(phaseIdx);
 
             const auto& tij = fluxVarsCache.advectionTij();
-            const auto& g = problem.gravityAtPos(scvf.ipGlobal());
+            const auto& g = problem.spatialParams().gravity(scvf.ipGlobal());
 
             //! compute alpha := n^T*K*g
             const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
@@ -350,7 +350,7 @@ public:
             }();
 
             const auto& tij = fluxVarsCache.advectionTij();
-            const auto& g = problem.gravityAtPos(scvf.ipGlobal());
+            const auto& g = problem.spatialParams().gravity(scvf.ipGlobal());
 
             // Obtain inside and outside pressures
             const auto pInside = insideVolVars.pressure(phaseIdx);
diff --git a/dumux/geomechanics/elastic/localresidual.hh b/dumux/geomechanics/elastic/localresidual.hh
index c0db87e23832da9b10980aef581215258e8c574d..dcad48af0b5db4f97165503c40f87b63b6113f75 100644
--- a/dumux/geomechanics/elastic/localresidual.hh
+++ b/dumux/geomechanics/elastic/localresidual.hh
@@ -128,7 +128,7 @@ public:
         static const bool gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
         if (gravity)
         {
-            const auto& g = problem.gravityAtPos(scv.center());
+            const auto& g = problem.spatialParams().gravity(scv.center());
             for (int dir = 0; dir < GridView::dimensionworld; ++dir)
                 source[Indices::momentum(dir)] += elemVolVars[scv].solidDensity()*g[dir];
         }
diff --git a/dumux/geomechanics/poroelastic/localresidual.hh b/dumux/geomechanics/poroelastic/localresidual.hh
index 387dd75fd53b8e39e049893fd89d7b414d2b8f85..14517966cd44971143393dafbd5f52703b5b5d0a 100644
--- a/dumux/geomechanics/poroelastic/localresidual.hh
+++ b/dumux/geomechanics/poroelastic/localresidual.hh
@@ -94,7 +94,7 @@ public:
             const auto rhoAverage = phi*rhoFluid + (1.0 - phi*vv.solidDensity());
 
             // add body force
-            const auto& g = problem.gravityAtPos(scv.center());
+            const auto& g = problem.spatialParams().gravity(scv.center());
             for (int dir = 0; dir < GridView::dimensionworld; ++dir)
                 source[ Indices::momentum(dir) ] += rhoAverage*g[dir];
         }
diff --git a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh
index 077f00c640ec2291c85f46444fa00402684962fd..de16ec1e6a545b35dffe06a5b5b7a3ca805ffb45 100644
--- a/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh
+++ b/dumux/multidomain/boundary/darcydarcy/couplingmanager.hh
@@ -186,7 +186,7 @@ public:
         {
             // do averaging for the density over all neighboring elements
             const auto rho = (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5;
-            const auto& g = this->problem(domainI).gravityAtPos(scvf.ipGlobal());
+            const auto& g = this->problem(domainI).spatialParams().gravity(scvf.ipGlobal());
 
             //! compute alpha := n^T*K*g
             const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor();
diff --git a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
index 854fa853bcdbf18d792115a7b59afc80db9db45c..47162178cd1de6d3671d92ab52fade7437228017 100644
--- a/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
+++ b/dumux/multidomain/boundary/stokesdarcy/couplingdata.hh
@@ -455,7 +455,7 @@ protected:
         const Scalar mu = context.volVars.viscosity(darcyPhaseIdx);
         const Scalar rho = context.volVars.density(darcyPhaseIdx);
         const Scalar distance = (context.element.geometry().center() - scvf.center()).two_norm();
-        const Scalar g = -scvf.directionSign() * couplingManager_.problem(darcyIdx).gravity()[scvf.directionIndex()];
+        const Scalar g = -scvf.directionSign() * couplingManager_.problem(darcyIdx).spatialParams().gravity(scvf.center())[scvf.directionIndex()];
         const Scalar interfacePressure = ((scvf.directionSign() * velocity * (mu/darcyPermeability(element, scvf))) + rho * g) * distance + cellCenterPressure;
         return interfacePressure;
     }
@@ -480,7 +480,7 @@ protected:
         const Scalar mu = context.volVars.viscosity(darcyPhaseIdx);
         const Scalar rho = context.volVars.density(darcyPhaseIdx);
         const Scalar distance = (context.element.geometry().center() - scvf.center()).two_norm();
-        const Scalar g = -scvf.directionSign() * couplingManager_.problem(darcyIdx).gravity()[scvf.directionIndex()];
+        const Scalar g = -scvf.directionSign() * couplingManager_.problem(darcyIdx).spatialParams().gravity(scvf.center())[scvf.directionIndex()];
 
         // get the Forchheimer coefficient
         Scalar cF = 0.0;
diff --git a/dumux/multidomain/facet/box/darcyslaw.hh b/dumux/multidomain/facet/box/darcyslaw.hh
index c2045337017cf82c35ead2bcda36f3fa3a781398..585bcd5106ea765db33f5d2474b6452f2202251d 100644
--- a/dumux/multidomain/facet/box/darcyslaw.hh
+++ b/dumux/multidomain/facet/box/darcyslaw.hh
@@ -123,7 +123,7 @@ public:
         static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity");
         if (enableGravity)
             flux -= rho*scvf.area()*insideVolVars.extrusionFactor()
-                    *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), problem.gravityAtPos(scvf.center()));
+                    *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), problem.spatialParams().gravity(scvf.center()));
 
         return flux;
     }
diff --git a/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh b/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh
index ac321eaf425c97286a6a503144662acb22852643..87d30154df7794de350b4ca03d2379ff2bd9a094 100644
--- a/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh
+++ b/dumux/multidomain/facet/cellcentered/mpfa/localassembler.hh
@@ -147,7 +147,7 @@ public:
             // gravitational acceleration on this face
             const auto& curLocalScvf = iv.localScvf(faceIdx);
             const auto& curGlobalScvf = this->fvGeometry().scvf(curLocalScvf.gridScvfIndex());
-            const auto& gravity = this->problem().gravityAtPos(curGlobalScvf.ipGlobal());
+            const auto& gravity = this->problem().spatialParams().gravity(curGlobalScvf.ipGlobal());
             const auto curIsInteriorBoundary = curLocalScvf.isOnInteriorBoundary();
             const Scalar curXiFactor = curIsInteriorBoundary ? (curGlobalScvf.boundary() ? 1.0 : xi) : 1.0;
 
diff --git a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh
index 2350e9fdbded6aa2ff5d0f412ac02ae169161ca4..fe782aba49f00d91c4bdd810a520826b936f3560 100644
--- a/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh
+++ b/dumux/multidomain/facet/cellcentered/tpfa/darcyslaw.hh
@@ -201,7 +201,7 @@ class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, FVGridGeometry, /*isNetwork*/
         if (gravity)
         {
             // compute alpha := n^T*K*g and add to flux (use arithmetic mean for density)
-            const auto& g = problem.gravityAtPos(scvf.ipGlobal());
+            const auto& g = problem.spatialParams().gravity(scvf.ipGlobal());
             const auto rho = 0.5*(insideVolVars.density(phaseIdx) + facetVolVars.density(phaseIdx));
             const auto rhoTimesArea = rho*scvf.area();
             const auto alpha_inside = rhoTimesArea*insideVolVars.extrusionFactor()
diff --git a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/problem_darcy.hh b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/problem_darcy.hh
index b7c966b640b69be3a891d27434d0d5df9af711c5..4da06f847293acccdf5065b1e77e9513dc73b841 100644
--- a/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/problem_darcy.hh
+++ b/test/multidomain/boundary/stokesdarcy/1p2c_2p2c/problem_darcy.hh
@@ -312,7 +312,7 @@ public:
         PrimaryVariables values(0.0);
         values.setState(initialPhasePresence_);
 
-        values[pressureIdx] = pressure_ + 1000. * this->gravity()[1] * (globalPos[1] - this->fvGridGeometry().bBoxMax()[1]);
+        values[pressureIdx] = pressure_ + 1000. * this->spatialParams().gravity(globalPos)[1] * (globalPos[1] - this->fvGridGeometry().bBoxMax()[1]);
         values[switchIdx] = initialSw_;
 
 #if NONISOTHERMAL
diff --git a/test/multidomain/facet/1p_1p/gravity/problem_bulk.hh b/test/multidomain/facet/1p_1p/gravity/problem_bulk.hh
index a029617b91129ceec18f1eb51d9386112d401cc2..c6efebf998f474a83500d45319eecf70837cdec2 100644
--- a/test/multidomain/facet/1p_1p/gravity/problem_bulk.hh
+++ b/test/multidomain/facet/1p_1p/gravity/problem_bulk.hh
@@ -162,7 +162,7 @@ public:
     //! Evaluates the initial conditions.
     PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
     {
-        const auto g = this->gravityAtPos(globalPos)[dimWorld-1];
+        const auto g = this->spatialParams().gravity(globalPos)[dimWorld-1];
         const auto h = 1.0 - (3.0-globalPos[dimWorld-1])*g;
         return PrimaryVariables(h);
     }
diff --git a/test/multidomain/facet/1p_1p/gravity/problem_lowdim.hh b/test/multidomain/facet/1p_1p/gravity/problem_lowdim.hh
index a2f1c82933bcfb1e1e3ffd43c620420f935ad348..fe61484a889b7c8d351981d70ba31431cdbd1b3e 100644
--- a/test/multidomain/facet/1p_1p/gravity/problem_lowdim.hh
+++ b/test/multidomain/facet/1p_1p/gravity/problem_lowdim.hh
@@ -158,7 +158,7 @@ public:
     //! Evaluates the initial conditions.
     PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
     {
-        const auto g = this->gravityAtPos(globalPos)[dimWorld-1];
+        const auto g = this->spatialParams().gravity(globalPos)[dimWorld-1];
         const auto h = 1.0 - (3.0-globalPos[dimWorld-1])*g;
         return PrimaryVariables(h);
     }
diff --git a/test/porousmediumflow/2p/implicit/cornerpoint/problem.hh b/test/porousmediumflow/2p/implicit/cornerpoint/problem.hh
index a11788dee5e3923ab4566bf18dacd7eb6b8400c4..89b4e89442640ef64362b683b9a3acc1abde4158 100644
--- a/test/porousmediumflow/2p/implicit/cornerpoint/problem.hh
+++ b/test/porousmediumflow/2p/implicit/cornerpoint/problem.hh
@@ -210,7 +210,7 @@ public:
 
         // hydrostatic pressure
         Scalar densityW = 1000;
-        values[Indices::pressureIdx] = 1e5 + densityW*(this->gravity()*globalPos);
+        values[Indices::pressureIdx] = 1e5 + densityW*(this->spatialParams().gravity(globalPos)*globalPos);
         values[Indices::saturationIdx] = 0.0;
 
         return values;
diff --git a/test/porousmediumflow/2p/implicit/fracture/problem.hh b/test/porousmediumflow/2p/implicit/fracture/problem.hh
index af92f08a24041532b3dc21beb919b1b9922abf73..3292340e70cefee3ceae2bac806c9116be5764a1 100644
--- a/test/porousmediumflow/2p/implicit/fracture/problem.hh
+++ b/test/porousmediumflow/2p/implicit/fracture/problem.hh
@@ -194,7 +194,7 @@ public:
     PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
     {
         const auto depth = this->fvGridGeometry().bBoxMax()[dimWorld-1] - globalPos[dimWorld-1];
-        const auto g = this->gravityAtPos(globalPos)[dimWorld-1];
+        const auto g = this->spatialParams().gravity(globalPos)[dimWorld-1];
 
         PrimaryVariables values;
         values[pressureIdx] = 1e5 + 1000*g*depth;
diff --git a/test/porousmediumflow/2p/implicit/incompressible/problem.hh b/test/porousmediumflow/2p/implicit/incompressible/problem.hh
index 2f99bf47b99f9e8ed06352008ba029e85838ef84..100b86651ede0a06703dd9917b32b058362c1984 100644
--- a/test/porousmediumflow/2p/implicit/incompressible/problem.hh
+++ b/test/porousmediumflow/2p/implicit/incompressible/problem.hh
@@ -171,7 +171,7 @@ public:
         Scalar factor = (width*alpha + (1.0 - alpha)*globalPos[0])/width;
 
         // hydrostatic pressure scaled by alpha
-        values[pressureH2OIdx] = 1e5 - factor*densityW*this->gravity()[1]*depth;
+        values[pressureH2OIdx] = 1e5 - factor*densityW*this->spatialParams().gravity(globalPos)[1]*depth;
         values[saturationDNAPLIdx] = 0.0;
 
         return values;
@@ -216,7 +216,7 @@ public:
         Scalar depth = this->fvGridGeometry().bBoxMax()[1] - globalPos[1];
 
         // hydrostatic pressure
-        values[pressureH2OIdx] = 1e5 - densityW*this->gravity()[1]*depth;
+        values[pressureH2OIdx] = 1e5 - densityW*this->spatialParams().gravity(globalPos)[1]*depth;
         values[saturationDNAPLIdx] = 0;
         return values;
     }
diff --git a/test/porousmediumflow/2p2c/implicit/injection/problem.hh b/test/porousmediumflow/2p2c/implicit/injection/problem.hh
index fe1c10d0ecb7a9738176621cb483b7afa5726a75..0a2073f4befa9a57a6e7f4c0e7a8f8652db41792 100644
--- a/test/porousmediumflow/2p2c/implicit/injection/problem.hh
+++ b/test/porousmediumflow/2p2c/implicit/injection/problem.hh
@@ -333,7 +333,7 @@ private:
 
         Scalar densityW = FluidSystem::H2O::liquidDensity(temperature_, 1e5);
 
-        Scalar pl = 1e5 - densityW*this->gravity()[1]*(depthBOR_ - globalPos[1]);
+        Scalar pl = 1e5 - densityW*this->spatialParams().gravity(globalPos)[1]*(depthBOR_ - globalPos[1]);
         Scalar moleFracLiquidN2 = pl*0.95/BinaryCoeff::H2O_N2::henry(temperature_);
         Scalar moleFracLiquidH2O = 1.0 - moleFracLiquidN2;
 
diff --git a/test/porousmediumflow/co2/implicit/problem.hh b/test/porousmediumflow/co2/implicit/problem.hh
index be01eb231c66129f4a5c40a8e877fa1ff58f52de..654652729447c062d65e9a0805fdd471000c67c6 100644
--- a/test/porousmediumflow/co2/implicit/problem.hh
+++ b/test/porousmediumflow/co2/implicit/problem.hh
@@ -474,7 +474,7 @@ private:
         else // mass-fraction formulation
             values[switchIdx] = moleFracLiquidCO2*FluidSystem::molarMass(CO2Idx)/meanM;
 
-        values[pressureIdx] = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(depthBOR_ - globalPos[dimWorld-1]);
+        values[pressureIdx] = 1.0e5 - densityW*this->spatialParams().gravity(globalPos)[dimWorld-1]*(depthBOR_ - globalPos[dimWorld-1]);
 
 #if !ISOTHERMAL
         values[temperatureIdx] = temp;