diff --git a/dumux/common/boundingboxtree.hh b/dumux/common/boundingboxtree.hh
deleted file mode 100644
index 9c6ed9623bc94576eab564c2b6cfe7b60f22328b..0000000000000000000000000000000000000000
--- a/dumux/common/boundingboxtree.hh
+++ /dev/null
@@ -1,21 +0,0 @@
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program 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.                                     *
- *                                                                           *
- *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-#ifndef DUMUX_OLD_BOUNDINGBOXTREE_HH
-#define DUMUX_OLD_BOUNDINGBOXTREE_HH
-#warning "This header is deprecated. Include dumux/common/geometry/boundingboxtree.hh"
-#include <dumux/common/geometry/boundingboxtree.hh>
-#endif
diff --git a/dumux/discretization/fvgridvariables.hh b/dumux/discretization/fvgridvariables.hh
index a5f05550cd05215e403652de32bfe1d2b3522468..ce693bd28360c4e92c0f8276d6485fde9fc2b2b0 100644
--- a/dumux/discretization/fvgridvariables.hh
+++ b/dumux/discretization/fvgridvariables.hh
@@ -85,15 +85,6 @@ public:
         prevGridVolVars_ = curGridVolVars_;
     }
 
-    //! initialize all variables (instationary case)
-    template<class SolutionVector>
-    [[deprecated("Use init(sol) instead. The class now works without modification for stationary and instationary cases.")]]
-    void init(const SolutionVector& curSol, const SolutionVector& initSol)
-    {
-        // initialize current volvars and the flux var cache
-        init(initSol);
-    }
-
     //! update all variables
     template<class SolutionVector>
     void update(const SolutionVector& curSol, bool forceFluxCacheUpdate = false)
diff --git a/dumux/discretization/staggered/fvgridgeometry.hh b/dumux/discretization/staggered/fvgridgeometry.hh
index 466b48f1ce122209028643a797e7414b342055ad..66faf2a07af3e80c3944e7bc1e1abb31c4ab8f47 100644
--- a/dumux/discretization/staggered/fvgridgeometry.hh
+++ b/dumux/discretization/staggered/fvgridgeometry.hh
@@ -96,13 +96,6 @@ public:
     const auto& elementMapper() const
     { return gridGeometry_->elementMapper(); }
 
-    /*!
-     * \brief Returns the actual gridGeometry we are a restriction of
-     */
-    [[deprecated("Use actualGridGeometry instead")]]
-    const ActualGridGeometry& actualfvGridGeometry() const
-    { return actualGridGeometry(); }
-
     /*!
      * \brief Returns the actual gridGeometry we are a restriction of
      */
diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh
index d40c1485a635ae6acb625b73b8763ac25277d4f4..4f495a88cd466466c8c25ae59a549554086bbe66 100644
--- a/dumux/multidomain/embedded/couplingmanager1d3d.hh
+++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh
@@ -127,15 +127,15 @@ public:
         for (const auto& is : intersections(this->glue()))
         {
             // all inside elements are identical...
-            const auto& inside = is.inside(0);
+            const auto& inside = is.targetEntity(0);
             const auto intersectionGeometry = is.geometry();
             const std::size_t lowDimElementIdx = this->problem(lowDimIdx).gridGeometry().elementMapper().index(inside);
 
             // compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
             const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx);
-            for (std::size_t outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
+            for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
             {
-                const auto& outside = is.outside(outsideIdx);
+                const auto& outside = is.domainEntity(outsideIdx);
                 const std::size_t bulkElementIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(outside);
                 lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
             }
@@ -500,15 +500,15 @@ public:
         for (const auto& is : intersections(this->glue()))
         {
             // all inside elements are identical...
-            const auto& inside = is.inside(0);
+            const auto& inside = is.targetEntity(0);
             const auto intersectionGeometry = is.geometry();
             const std::size_t lowDimElementIdx = this->problem(lowDimIdx).gridGeometry().elementMapper().index(inside);
 
             // compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
             const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx);
-            for (std::size_t outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
+            for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
             {
-                const auto& outside = is.outside(outsideIdx);
+                const auto& outside = is.domainEntity(outsideIdx);
                 const std::size_t bulkElementIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(outside);
                 lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
             }
@@ -778,15 +778,15 @@ public:
         for (const auto& is : intersections(this->glue()))
         {
             // all inside elements are identical...
-            const auto& inside = is.inside(0);
+            const auto& inside = is.targetEntity(0);
             const auto intersectionGeometry = is.geometry();
             const std::size_t lowDimElementIdx = this->problem(lowDimIdx).gridGeometry().elementMapper().index(inside);
 
             // compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
             const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx);
-            for (std::size_t outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
+            for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
             {
-                const auto& outside = is.outside(outsideIdx);
+                const auto& outside = is.domainEntity(outsideIdx);
                 const std::size_t bulkElementIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(outside);
                 lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
             }
@@ -952,7 +952,7 @@ public:
         for (const auto& is : intersections(this->glue()))
         {
             // all inside elements are identical...
-            const auto& inside = is.inside(0);
+            const auto& inside = is.targetEntity(0);
             // get the intersection geometry for integrating over it
             const auto intersectionGeometry = is.geometry();
 
@@ -966,9 +966,9 @@ public:
             for (auto&& qp : quad)
             {
                 // compute the coupling stencils
-                for (std::size_t outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
+                for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
                 {
-                    const auto& outside = is.outside(outsideIdx);
+                    const auto& outside = is.domainEntity(outsideIdx);
                     const std::size_t bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside);
 
                     // each quadrature point will be a point source for the sub problem
@@ -977,8 +977,8 @@ public:
                     const auto qpweight = qp.weight();
                     const auto ie = intersectionGeometry.integrationElement(qp.position());
                     this->pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, ie, std::vector<std::size_t>({lowDimElementIdx}));
-                    this->pointSources(lowDimIdx).back().setEmbeddings(is.neighbor(0));
-                    computeBulkSource(globalPos, kernelWidth, id, lowDimElementIdx, bulkElementIdx, qpweight*ie/is.neighbor(0));
+                    this->pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors());
+                    computeBulkSource(globalPos, kernelWidth, id, lowDimElementIdx, bulkElementIdx, qpweight*ie/is.numDomainNeighbors());
 
                     // pre compute additional data used for the evaluation of
                     // the actual solution dependent source term
@@ -1083,15 +1083,15 @@ public:
         for (const auto& is : intersections(this->glue()))
         {
             // all inside elements are identical...
-            const auto& inside = is.inside(0);
+            const auto& inside = is.targetEntity(0);
             const auto intersectionGeometry = is.geometry();
             const std::size_t lowDimElementIdx = this->problem(lowDimIdx).gridGeometry().elementMapper().index(inside);
 
             // compute the volume the low-dim domain occupies in the bulk domain if it were full-dimensional
             const auto radius = this->problem(lowDimIdx).spatialParams().radius(lowDimElementIdx);
-            for (std::size_t outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
+            for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
             {
-                const auto& outside = is.outside(outsideIdx);
+                const auto& outside = is.domainEntity(outsideIdx);
                 const std::size_t bulkElementIdx = this->problem(bulkIdx).gridGeometry().elementMapper().index(outside);
                 lowDimVolumeInBulkElement_[bulkElementIdx] += intersectionGeometry.volume()*M_PI*radius*radius;
             }
diff --git a/dumux/multidomain/embedded/couplingmanagerbase.hh b/dumux/multidomain/embedded/couplingmanagerbase.hh
index f64b47db970ded62e2d280f56840758750886c00..678eb0ad4ce9b69dbd97a6d3f1475a0564f1f15b 100644
--- a/dumux/multidomain/embedded/couplingmanagerbase.hh
+++ b/dumux/multidomain/embedded/couplingmanagerbase.hh
@@ -257,7 +257,7 @@ public:
         for (const auto& is : intersections(this->glue()))
         {
             // all inside elements are identical...
-            const auto& inside = is.inside(0);
+            const auto& inside = is.targetEntity(0);
             // get the intersection geometry for integrating over it
             const auto intersectionGeometry = is.geometry();
 
@@ -269,9 +269,9 @@ public:
             for (auto&& qp : quad)
             {
                 // compute the coupling stencils
-                for (std::size_t outsideIdx = 0; outsideIdx < is.neighbor(0); ++outsideIdx)
+                for (std::size_t outsideIdx = 0; outsideIdx < is.numDomainNeighbors(); ++outsideIdx)
                 {
-                    const auto& outside = is.outside(outsideIdx);
+                    const auto& outside = is.domainEntity(outsideIdx);
                     const std::size_t bulkElementIdx = bulkFvGridGeometry.elementMapper().index(outside);
 
                     // each quadrature point will be a point source for the sub problem
@@ -280,9 +280,9 @@ public:
                     const auto qpweight = qp.weight();
                     const auto ie = intersectionGeometry.integrationElement(qp.position());
                     pointSources(bulkIdx).emplace_back(globalPos, id, qpweight, ie, std::vector<std::size_t>({bulkElementIdx}));
-                    pointSources(bulkIdx).back().setEmbeddings(is.neighbor(0));
+                    pointSources(bulkIdx).back().setEmbeddings(is.numDomainNeighbors());
                     pointSources(lowDimIdx).emplace_back(globalPos, id, qpweight, ie, std::vector<std::size_t>({lowDimElementIdx}));
-                    pointSources(lowDimIdx).back().setEmbeddings(is.neighbor(0));
+                    pointSources(lowDimIdx).back().setEmbeddings(is.numDomainNeighbors());
 
                     // pre compute additional data used for the evaluation of
                     // the actual solution dependent source term
diff --git a/dumux/multidomain/glue.hh b/dumux/multidomain/glue.hh
index a8af2b8b62ba7be31e2c2c39d308592a357b385e..dab4af6316f7e194ac471491928bd58642d8d835 100644
--- a/dumux/multidomain/glue.hh
+++ b/dumux/multidomain/glue.hh
@@ -153,21 +153,6 @@ public:
     TargetElement targetEntity(unsigned int n = 0) const
     { return targetTree_.entitySet().entity(std::get<targetIdx>(neighbors_)[n]); }
 
-    //! get the number of neigbors
-    [[deprecated("neighbor is deprecated and will be removed after release 3.1. Use numDomainNeighbors() / numTargetNeighbors()")]]
-    std::size_t neighbor(unsigned int side) const
-    { return std::max(std::get<domainIdx>(neighbors_).size(), std::get<targetIdx>(neighbors_).size()); }
-
-    //! get the inside (domain) neighbor
-    [[deprecated("outside is deprecated and will be removed after release 3.1. Use domainIdx(uint)")]]
-    DomainElement outside(unsigned int n) const
-    { return domainTree_.entitySet().entity(std::get<domainIdx>(neighbors_)[n]); }
-
-    //! get the outside (target) neighbor
-    [[deprecated("inside is deprecated and will be removed after release 3.1. Use targetIdx(uint)")]]
-    TargetElement inside(unsigned int n) const
-    { return targetTree_.entitySet().entity(std::get<targetIdx>(neighbors_)[n]); }
-
 private:
     IndexStorage neighbors_;
     std::vector<GlobalPosition> corners_;
diff --git a/dumux/porousmediumflow/2p2c/sequential/celldata.hh b/dumux/porousmediumflow/2p2c/sequential/celldata.hh
index d26ed2fbec549fce031791bc5382aa6b8cd04555..b9ed61faa25642eddbdb5e602442902b35aa4c0d 100644
--- a/dumux/porousmediumflow/2p2c/sequential/celldata.hh
+++ b/dumux/porousmediumflow/2p2c/sequential/celldata.hh
@@ -128,17 +128,6 @@ public:
         return fluidState_->pressure(phaseIdx);
     }
 
-    /*!
-     * \brief Modify the phase pressure
-     * \param phaseIdx index of the Phase
-     * \param value Value to be stored
-     */
-    [[deprecated("cellData.setPressure is deprecated and will be removed after release 3.1. Use cellData.manipulateFluidState() and set pressure in fluid state directly")]]
-    void setPressure(int phaseIdx, Scalar value)
-    {
-        manipulateFluidState().setPressure(phaseIdx, value);
-    }
-
     /*!
      * \brief Returns the total mass concentration of a component \f$\mathrm{[kg/m^3]}\f$.
      *
@@ -294,26 +283,12 @@ public:
 
     /*** b) from fluidstate ***/
 
-    //! DOC ME!
-    [[deprecated("cellData.setSaturation is deprecated and will be removed after release 3.1. Use cellData.manipulateFluidState() and set saturation in fluid state directly")]]
-    void setSaturation(int phaseIdx, Scalar value)
-    {
-        fluidState_->setSaturation(phaseIdx, value);
-        fluidState_->setSaturation(1-phaseIdx, 1.0-value);
-    }
-
     //! DOC ME!
     const Scalar saturation(int phaseIdx) const
     {
         return fluidState_->saturation(phaseIdx);
     }
 
-    //! DOC ME!
-    [[deprecated("cellData.setViscosity is deprecated and will be removed after release 3.1. Use cellData.manipulateFluidState() and set viscosity in fluid state directly")]]
-    void setViscosity(int phaseIdx, Scalar value)
-    {
-        fluidState_->setViscosity(phaseIdx, value);
-    }
     //! DOC ME!
     const Scalar viscosity(int phaseIdx) const
     {