diff --git a/CHANGELOG.md b/CHANGELOG.md
index c3f3c89a428db01d7ea95077f4689fc1afa85f28..39f99338b5413b7e4aa39d03f54c06dbc0bc50fb 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -32,6 +32,7 @@ idenpendently of the property system.
 - __Thermal conductivity__: Removed deprecated alias `ThermalConductivitySomerton`. Use `ThermalConductivitySomertonTwoP` and `ThermalConductivitySomertonThreeP` instead.
 - __Thermal conductivity__: `1p` folder from `material/fluidmatrixinteractions` that was deprecated in 3.9 has been removed. Use `ThermalConductivityAverage` from `dumux/material/fluidmatrixinteractions/thermalconductivityaverage.hh` instead.
 - __setCheckPoint(vec)__: this function in `timeloop.hh` is removed. Other options are available.
+- __BJS functions__: the functions in `freeflow/navierstokes/momentum/problem.hh` returning permeability, alpha, beta, and pm velocity to be used in BJS boundary condition are removed. They are supposed to be implemented in the test problem.
 
 Differences Between DuMu<sup>x</sup> 3.9 and DuMu<sup>x</sup> 3.8
 =============================================
diff --git a/dumux/freeflow/navierstokes/momentum/problem.hh b/dumux/freeflow/navierstokes/momentum/problem.hh
index 4f560a277a8bc00158938f086dd5621d6b1174ad..37bcbc51147d340c05ec5724e3f8c690457b6b14 100644
--- a/dumux/freeflow/navierstokes/momentum/problem.hh
+++ b/dumux/freeflow/navierstokes/momentum/problem.hh
@@ -423,89 +423,6 @@ public:
     bool onSlipBoundaryAtPos(const GlobalPosition& pos) const
     { return false; }
 
-    /*!
-     * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
-     * This member function must be overloaded in the problem implementation, if the BJS boundary condition is used.
-     */
-    [[deprecated("Will be removed after release 3.9.")]]
-    Scalar permeability(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
-    {
-        DUNE_THROW(Dune::NotImplemented,
-            "When using the Beavers-Joseph-Saffman boundary condition, "
-            "the permeability must be returned in the actual problem"
-        );
-    }
-
-    /*!
-     * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
-     * This member function must be overloaded in the problem implementation, if the BJS boundary condition is used.
-     */
-    [[deprecated("Will be removed after release 3.9. Implement betaBJ instead. ")]]
-    Scalar alphaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
-    {
-        DUNE_THROW(Dune::NotImplemented,
-            "When using the Beavers-Joseph-Saffman boundary condition, "
-            "the alpha value must be returned in the actual problem"
-        );
-    }
-
-    /*!
-     * \brief Returns the beta value which is the alpha value divided by the square root of the (scalar-valued) interface permeability.
-     */
-    [[deprecated("Needs to be implemented in test problem. Will be removed after release 3.9.")]]
-    Scalar betaBJ(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf, const GlobalPosition& tangentialVector) const
-    {
-        const auto& K = asImp_().permeability(fvGeometry, scvf);
-        const auto interfacePermeability = vtmv(tangentialVector, K, tangentialVector);
-        using std::sqrt;
-        return asImp_().alphaBJ(fvGeometry, scvf) / sqrt(interfacePermeability);
-    }
-
-    /*!
-     * \brief Returns the velocity in the porous medium (which is 0 by default according to Saffmann).
-     */
-    [[deprecated("Needs to be implemented in test problem. Will be removed after release 3.9.")]]
-    VelocityVector porousMediumVelocity(const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) const
-    {
-        return VelocityVector(0.0);
-    }
-
-    /*!
-     * \brief Returns the slip velocity at a porous boundary based on the Beavers-Joseph(-Saffman) condition.
-     * \note This only returns a vector filled with one component of the slip velocity (corresponding to the dof axis of the scv the scvf belongs to)
-     */
-    [[deprecated("Will be removed after release 3.9.")]]
-    const VelocityVector beaversJosephVelocity(const FVElementGeometry& fvGeometry,
-                                               const SubControlVolumeFace& scvf,
-                                               const ElementVolumeVariables& elemVolVars,
-                                               Scalar tangentialVelocityGradient) const
-    {
-        assert(scvf.isLateral());
-        assert(scvf.boundary());
-
-        const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
-
-        // create a unit normal vector oriented in positive coordinate direction
-        GlobalPosition orientation(0.0);
-        orientation[scv.dofAxis()] = 1.0;
-
-        // du/dy + dv/dx = beta * (u_boundary-uPM)
-        // beta = alpha/sqrt(K)
-        const Scalar betaBJ = asImp_().betaBJ(fvGeometry, scvf, orientation);
-        const Scalar distanceNormalToBoundary = (scvf.ipGlobal() - scv.dofPosition()).two_norm();
-
-        static const bool onlyNormalGradient = getParamFromGroup<bool>(this->paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradientForBeaversJoseph", false);
-        if (onlyNormalGradient)
-            tangentialVelocityGradient = 0.0;
-
-        const Scalar scalarSlipVelocity = (tangentialVelocityGradient*distanceNormalToBoundary
-            + asImp_().porousMediumVelocity(fvGeometry, scvf) * orientation * betaBJ * distanceNormalToBoundary
-            + elemVolVars[scv].velocity()) / (betaBJ*distanceNormalToBoundary + 1.0);
-
-        orientation *= scalarSlipVelocity;
-        return orientation;
-    }
-
     const CouplingManager& couplingManager() const
     {
         if constexpr (isCoupled_)