diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh
index 292ddf3be65c2ca6438e3881c1e499af0071e15b..7cec280faf7aeeb7f5ad4be2c359031e75f47912 100644
--- a/dumux/assembly/boxlocalassembler.hh
+++ b/dumux/assembly/boxlocalassembler.hh
@@ -560,7 +560,7 @@ public:
                                                        problem,
                                                        element,
                                                        fvGeometry,
-                                                       volVars,
+                                                       curElemVolVars,
                                                        scv);
         }
 
diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh
index 9a9b8168e246e8cc6ae78c09022695a96db0926d..6855b933e86a0a75603753ea69160ec763be70c8 100644
--- a/dumux/assembly/cclocalassembler.hh
+++ b/dumux/assembly/cclocalassembler.hh
@@ -540,7 +540,7 @@ public:
             this->localResidual().addStorageDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
 
         // add source term derivatives
-        this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
+        this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, curElemVolVars, scv);
 
         // add flux derivatives for each scvf
         for (const auto& scvf : scvfs(fvGeometry))
diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index f3063569475d9ea073ead553445ac21c90927da7..49118f70bfe88bc82a554f21fb6f10c94bb771c3 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -415,11 +415,11 @@ public:
      * \brief Add source term derivative to the Jacobian
      * \note Only needed in case of analytic differentiation and solution dependent sources
      */
-    template<class MatrixBlock, class VolumeVariables>
+    template<class MatrixBlock, class ElementVolumeVariables>
     void addSourceDerivatives(MatrixBlock& block,
                               const Element& element,
                               const FVElementGeometry& fvGeometry,
-                              const VolumeVariables& volVars,
+                              const ElementVolumeVariables& curElemVolVars,
                               const SubControlVolume& scv) const {}
 
     /*!
diff --git a/dumux/multidomain/couplingmanager.hh b/dumux/multidomain/couplingmanager.hh
index b1fbbb629653c319ee4c3ab836b4586e279691d8..4583039bc0223eed06766dc4b698f9e4c6a82cc5 100644
--- a/dumux/multidomain/couplingmanager.hh
+++ b/dumux/multidomain/couplingmanager.hh
@@ -239,6 +239,28 @@ public:
         return NumericEpsilon<typename Traits::Scalar, numEq>(paramGroup);
     }
 
+    /*!
+     * \ingroup MultiDomain
+     * \brief Evaluates the element residual derivative of a coupled element of domain i
+     *        with respect to all coupled degrees of freedom of domain j (for analytic Jacobians)
+     *
+     * \param Aij the coupling block of Jacobian between domain i and j
+     * \param domainI the domain index of domain i
+     * \param localAssemblerI the local assembler assembling the element residual of an element of domain i
+     * \param domainJ the domain index of domain j
+     *
+     * \note  The element whose residual derivatives are to be evaluated can be retrieved from the local assembler
+     *        as localAssemblerI.element() as well as all up-to-date variables and caches.
+     */
+    template<std::size_t i, std::size_t j, class LocalAssemblerI, class JacobianCouplingBlock>
+    void addCouplingDerivatives(JacobianCouplingBlock& Aij,
+                                Dune::index_constant<i> domainI,
+                                const LocalAssemblerI& localAssemblerI,
+                                Dune::index_constant<j> domainJ) const
+    {
+        DUNE_THROW(Dune::NotImplemented, "This coupling manager does not implement analytic derivatives. Use numeric differentiation instead.");
+    }
+
     /*!
      * \brief set the pointers to the sub problems
      * \param problems A tuple of shared pointers to the sub problems
diff --git a/dumux/multidomain/subdomainboxlocalassembler.hh b/dumux/multidomain/subdomainboxlocalassembler.hh
index e70259f712fa839bda287c2923c1d11470e8e1cf..3b599e8e032391065ab9adee8b87c9a114d3cd7e 100644
--- a/dumux/multidomain/subdomainboxlocalassembler.hh
+++ b/dumux/multidomain/subdomainboxlocalassembler.hh
@@ -693,6 +693,296 @@ public:
     {}
 };
 
+/*!
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief Box local assembler using analytic differentiation and implicit time discretization
+ */
+template<std::size_t id, class TypeTag, class Assembler>
+class SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true>
+: public SubDomainBoxLocalAssemblerBase<id, TypeTag, Assembler,
+                                        SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, true>, true>
+{
+    using ThisType = SubDomainBoxLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, true>;
+    using ParentType = SubDomainBoxLocalAssemblerBase<id, TypeTag, Assembler, ThisType, true>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+    using ElementResidualVector = typename GetPropType<TypeTag, Properties::LocalResidual>::ElementResidualVector;
+    static constexpr auto domainI = Dune::index_constant<id>();
+    static constexpr int numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
+    static constexpr bool enableGridFluxVarsCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
+    static constexpr bool enableGridVolVarsCache = getPropValue<TypeTag, Properties::EnableGridVolumeVariablesCache>();
+    enum { dim = GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension };
+
+public:
+    using ParentType::ParentType;
+
+    /*!
+     * \brief Computes the derivatives with respect to the given element and adds them
+     *        to the global matrix.
+     *
+     * \return The element residual at the current solution.
+     */
+    ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
+    {
+        // get some aliases for convenience
+        const auto& element = this->element();
+        const auto& fvGeometry = this->fvGeometry();
+        const auto& curSol = this->curSol()[domainI];
+        const auto& problem = this->problem();
+        auto&& curElemVolVars = this->curElemVolVars();
+        const auto& elemFluxVarsCache = this->elemFluxVarsCache();
+
+        // get the vecor of the acutal element residuals
+        const auto origResiduals = this->evalLocalResidual();
+        auto origSourceResidual = origResiduals; origSourceResidual = 0.0;
+        static const bool useNumDiff = getParam<bool>("Assembly.UseNumericDifferentiationForSourceDerivatives", false);
+        if (useNumDiff)
+            origSourceResidual = this->evalLocalSourceResidual(element, curElemVolVars);
+
+        //////////////////////////////////////////////////////////////////////////////////////////////////
+        //                                                                                              //
+        // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
+        // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
+        // actual element. In the actual element we evaluate the derivative of the entire residual.     //
+        //                                                                                              //
+        //////////////////////////////////////////////////////////////////////////////////////////////////
+        // create the element solution
+        auto elemSol = elementSolution(element, curSol, fvGeometry.fvGridGeometry());
+
+        // create the vector storing the partial derivatives
+        auto partialDerivs = origResiduals;
+        partialDerivs = 0.0;
+
+        // calculation of the source and storage derivatives
+        for (const auto& scv : scvs(fvGeometry))
+        {
+            // dof index and corresponding actual pri vars
+            const auto dofIdx = scv.dofIndex();
+            const auto& volVars = curElemVolVars[scv];
+
+            // derivative of this scv residual w.r.t the d.o.f. of the same scv (because of mass lumping)
+            // only if the problem is instationary we add derivative of storage term
+            // TODO if e.g. porosity depends on all dofs in the element, we would have off-diagonal matrix entries!?
+            if (!this->assembler().isStationaryProblem())
+                this->localResidual().addStorageDerivatives(A[dofIdx][dofIdx], problem, element, fvGeometry, volVars, scv);
+
+            // derivative of this scv residual w.r.t the d.o.f. of the same scv (because of mass lumping)
+            // add source term derivatives
+            if (useNumDiff)
+            {
+                auto& curVolVars = this->getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
+                const auto origVolVars = curVolVars;
+
+                // calculate derivatives w.r.t to the privars at the dof at hand
+                for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+                {
+                    partialDerivs = 0.0;
+
+                    auto evalSource = [&](Scalar priVar)
+                    {
+                        const auto localDofIndex = scv.localDofIndex();
+                        elemSol[localDofIndex][pvIdx] = priVar;
+                        curVolVars.update(elemSol, this->problem(), element, scv);
+                        this->couplingManager().updateCouplingContext(domainI, *this, domainI, scv.dofIndex(), elemSol[localDofIndex], pvIdx);
+                        return this->evalLocalSourceResidual(element, curElemVolVars);
+                    };
+
+                    // derive the residuals numerically
+                    static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
+                    static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod");
+                    NumericDifferentiation::partialDerivative(evalSource, elemSol[scv.localDofIndex()][pvIdx], partialDerivs, origSourceResidual,
+                                                              eps_(elemSol[scv.localDofIndex()][pvIdx], pvIdx), numDiffMethod);
+
+                    // update the global stiffness matrix with the current partial derivatives
+                    for (auto&& scvJ : scvs(fvGeometry))
+                    {
+                        for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+                        {
+                            // A[i][col][eqIdx][pvIdx] is the rate of change of
+                            // the residual of equation 'eqIdx' at dof 'i'
+                            // depending on the primary variable 'pvIdx' at dof
+                            // 'col'.
+                            A[scvJ.dofIndex()][dofIdx][eqIdx][pvIdx] += partialDerivs[scvJ.localDofIndex()][eqIdx];
+                        }
+                    }
+
+                    // restore the original state of the scv's volume variables
+                    curVolVars = origVolVars;
+
+                    // restore the original element solution
+                    elemSol[scv.localDofIndex()][pvIdx] = curSol[scv.dofIndex()][pvIdx];
+
+                    // restore coupling manager
+                    this->couplingManager().updateCouplingContext(domainI, *this, domainI, scv.dofIndex(), elemSol[scv.localDofIndex()], pvIdx);
+                }
+            }
+            else
+                this->localResidual().addSourceDerivatives(A[dofIdx][dofIdx], problem, element, fvGeometry, curElemVolVars, scv);
+        }
+
+        // localJacobian[scvIdx][otherScvIdx][eqIdx][priVarIdx] of the fluxes
+        for (const auto& scvf : scvfs(fvGeometry))
+        {
+            // add flux term derivatives
+            if (!scvf.boundary())
+                this->localResidual().addFluxDerivatives(A, problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
+
+            // the boundary gets special treatment to simplify
+            // for the user
+            else
+            {
+                // add flux term derivatives
+                const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+                if (this->elemBcTypes()[insideScv.localDofIndex()].hasNeumann())
+                    this->localResidual().addRobinFluxDerivatives(A[insideScv.dofIndex()], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
+            }
+        }
+
+        // evaluate additional derivatives that might arise from the coupling
+        this->couplingManager().evalAdditionalDomainDerivatives(domainI, *this, origResiduals, A, gridVariables);
+
+        return origResiduals;
+    }
+
+    /*!
+     * \brief Computes the derivatives with respect to the given element and adds them
+     *        to the global matrix.
+     *
+     * \return The element residual at the current solution.
+     */
+    template<std::size_t otherId, class JacobianBlock, class GridVariables>
+    void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
+                                  const ElementResidualVector& res, GridVariables& gridVariables)
+    {
+        ///////////////////////////////////////////////////////
+        // Calculate derivatives of all dofs in the element  //
+        // with respect to all dofs in the coupling stencil. //
+        ///////////////////////////////////////////////////////
+        static const bool useNumDiff = getParam<bool>("Assembly.UseNumericDifferentiationForCouplingDerivatives", false);
+        if (useNumDiff)
+            assembleJacobianCouplingNumeric(domainJ, A, res, gridVariables);
+        else
+            this->couplingManager().addCouplingDerivatives(A, domainI, *this, domainJ);
+    }
+
+    /*!
+     * \brief Computes the derivatives with respect to the given element and adds them
+     *        to the global matrix.
+     * \todo This is just a copy of the assembleJacobianCoupling of the numeric-diff local assembler
+     * \return The element residual at the current solution.
+     */
+    template<std::size_t otherId, class JacobianBlock, class GridVariables>
+    void assembleJacobianCouplingNumeric(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
+                                         const ElementResidualVector& res, GridVariables& gridVariables)
+    {
+        ////////////////////////////////////////////////////////////////////////////////////////////////////////
+        // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
+        ////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+        // get some aliases for convenience
+        const auto& element = this->element();
+        const auto& fvGeometry = this->fvGeometry();
+        auto&& curElemVolVars = this->curElemVolVars();
+        auto&& elemFluxVarsCache = this->elemFluxVarsCache();
+
+        // get element stencil informations
+        const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ);
+
+        // convenience lambda for call to update self
+        auto updateSelf = [&] ()
+        {
+            // Update ourself after the context has been modified. Depending on the
+            // type of caching, other objects might have to be updated. All ifs can be optimized away.
+            if (enableGridFluxVarsCache)
+            {
+                if (enableGridVolVarsCache)
+                    this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache());
+                else
+                    this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, gridVariables.gridFluxVarsCache());
+            }
+            else
+            {
+                if (enableGridVolVarsCache)
+                    this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache);
+                else
+                    this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache);
+            }
+        };
+
+        const auto& curSolJ = this->curSol()[domainJ];
+        for (const auto globalJ : stencil)
+        {
+            // undeflected privars and privars to be deflected
+            const auto origPriVarsJ = curSolJ[globalJ];
+            auto priVarsJ = origPriVarsJ;
+
+            // the undeflected coupling residual
+            const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ);
+
+            for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx)
+            {
+                auto evalCouplingResidual = [&](Scalar priVar)
+                {
+                    priVarsJ[pvIdx] = priVar;
+                    this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
+                    updateSelf();
+                    return this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ);
+                };
+
+                // derive the residuals numerically
+                ElementResidualVector partialDerivs(element.subEntities(dim));
+
+                const auto& paramGroup = this->assembler().problem(domainJ).paramGroup();
+                static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod");
+                static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup);
+
+                NumericDifferentiation::partialDerivative(evalCouplingResidual, origPriVarsJ[pvIdx], partialDerivs, origResidual,
+                                                          epsCoupl(origPriVarsJ[pvIdx], pvIdx), numDiffMethod);
+
+                // update the global stiffness matrix with the current partial derivatives
+                for (auto&& scv : scvs(fvGeometry))
+                {
+                    for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+                    {
+                        // A[i][col][eqIdx][pvIdx] is the rate of change of
+                        // the residual of equation 'eqIdx' at dof 'i'
+                        // depending on the primary variable 'pvIdx' at dof
+                        // 'col'.
+
+                        // If the dof is coupled by a Dirichlet condition,
+                        // set the derived value only once (i.e. overwrite existing values).
+                        // For other dofs, add the contribution of the partial derivative.
+                        const auto bcTypes = this->elemBcTypes()[scv.localDofIndex()];
+                        if (bcTypes.isCouplingDirichlet(eqIdx))
+                            A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = partialDerivs[scv.localDofIndex()][eqIdx];
+                        else if (bcTypes.isDirichlet(eqIdx))
+                            A[scv.dofIndex()][globalJ][eqIdx][pvIdx] = 0.0;
+                        else
+                            A[scv.dofIndex()][globalJ][eqIdx][pvIdx] += partialDerivs[scv.localDofIndex()][eqIdx];
+                    }
+                }
+
+                // restore the current element solution
+                priVarsJ[pvIdx] = origPriVarsJ[pvIdx];
+
+                // restore the undeflected state of the coupling context
+                this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
+            }
+        }
+
+        // restore original state of the flux vars cache and/or vol vars in case of global caching.
+        // This has to be done in order to guarantee that everything is in an undeflected
+        // state before the assembly of another element is called. In the case of local caching
+        // this is obsolete because the local views used here go out of scope after this.
+        // We only have to do this for the last primary variable, for all others the flux var cache
+        // is updated with the correct element volume variables before residual evaluations
+        updateSelf();
+    }
+
+}; // implicit assebler with analytic Jacobian
+
 } // end namespace Dumux
 
 #endif
diff --git a/dumux/multidomain/subdomaincclocalassembler.hh b/dumux/multidomain/subdomaincclocalassembler.hh
index 12bd85cde509f826da298afd5c710cb876421fec..3882f22361d7ec7e60379bbf463b7fae9ba2da69 100644
--- a/dumux/multidomain/subdomaincclocalassembler.hh
+++ b/dumux/multidomain/subdomaincclocalassembler.hh
@@ -120,7 +120,7 @@ public:
 
         // for the diagonal jacobian block
         const auto globalI = this->fvGeometry().gridGeometry().elementMapper().index(this->element());
-        res[globalI] = this->asImp_().assembleJacobianAndResidualImplInverse(jacRow[domainId], *std::get<domainId>(gridVariables));
+        res[globalI] = this->asImp_().assembleJacobianAndResidualImpl(jacRow[domainId], *std::get<domainId>(gridVariables));
 
         // for the coupling blocks
         using namespace Dune::Hybrid;
@@ -308,7 +308,7 @@ public:
      * \return The element residual at the current solution.
      */
     template<class JacobianMatrixDiagBlock, class GridVariables>
-    LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
+    LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
     {
         //////////////////////////////////////////////////////////////////////////////////////////////////
         // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
@@ -640,7 +640,7 @@ public:
      * \return The element residual at the current solution.
      */
     template<class JacobianMatrixDiagBlock, class GridVariables>
-    LocalResidualValues assembleJacobianAndResidualImplInverse(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
+    LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
     {
         if (this->assembler().isStationaryProblem())
             DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
@@ -771,10 +771,12 @@ class SubDomainCCLocalAssembler<id, TypeTag, Assembler, DiffMethod::analytic, /*
     using Element = typename GridView::template Codim<0>::Entity;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
 
-    enum { numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq() };
-    enum { dim = GridView::dimension };
+    static constexpr int numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
+    static constexpr int dim = GridView::dimension;
 
     static constexpr auto domainI = Dune::index_constant<id>();
+    static constexpr bool enableGridFluxVarsCache = getPropValue<TypeTag, Properties::EnableGridFluxVariablesCache>();
+    static constexpr bool enableGridVolVarsCache = getPropValue<TypeTag, Properties::EnableGridVolumeVariablesCache>();
 
 public:
     using ParentType::ParentType;
@@ -792,11 +794,12 @@ public:
         const auto& problem = this->problem();
         const auto& element = this->element();
         const auto& fvGeometry = this->fvGeometry();
-        const auto& curElemVolVars = this->curElemVolVars();
+        const auto& gridGeometry = fvGeometry.gridGeometry();
+        auto&& curElemVolVars = this->curElemVolVars();
         const auto& elemFluxVarsCache = this->elemFluxVarsCache();
 
         // get reference to the element's current vol vars
-        const auto globalI = this->assembler().gridGeometry(domainI).elementMapper().index(element);
+        const auto globalI = gridGeometry.elementMapper().index(element);
         const auto& scv = fvGeometry.scv(globalI);
         const auto& volVars = curElemVolVars[scv];
 
@@ -805,7 +808,80 @@ public:
             this->localResidual().addStorageDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
 
         // add source term derivatives
-        this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
+        static const bool useNumDiff = getParam<bool>("Assembly.UseNumericDifferentiationForSourceDerivatives", false);
+        if (useNumDiff)
+        {
+            const auto sourceResidual = this->evalLocalSourceResidual(element, curElemVolVars)[0];
+            auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
+
+            // save a copy of the original privars and vol vars in order
+            // to restore the original solution after deflection
+            const auto& curSol = this->curSol()[domainI];
+            const auto origPriVars = curSol[globalI];
+            const auto origVolVars = curVolVars;
+
+            // element solution container to be deflected
+            auto elemSol = elementSolution(element, curSol, gridGeometry);
+
+            // derivatives in the neighbors with repect to the current elements
+            LocalResidualValues partialDeriv;
+            for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+            {
+                // reset derivatives of element dof with respect to itself
+                partialDeriv = 0.0;
+                auto evalSource = [&](Scalar priVar)
+                {
+                    // update the volume variables and calculate
+                    // the residual with the deflected primary variables
+                    elemSol[0][pvIdx] = priVar;
+                    this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalI, elemSol[0], pvIdx);
+                    curVolVars.update(elemSol, this->problem(), element, scv);
+                    return this->evalLocalSourceResidual(element, curElemVolVars)[0];
+                };
+
+                // for non-ghosts compute the derivative numerically
+                if (!this->elementIsGhost())
+                {
+                    static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
+                    static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod");
+                    NumericDifferentiation::partialDerivative(evalSource, elemSol[0][pvIdx], partialDeriv, sourceResidual,
+                                                              eps_(elemSol[0][pvIdx], pvIdx), numDiffMethod);
+                }
+
+                // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else
+                // as we always solve for a delta of the solution with repect to the initial solution this
+                // results in a delta of zero for ghosts
+                else partialDeriv[pvIdx] = 1.0;
+
+                // add the current partial derivatives to the global jacobian matrix
+                // only diagonal entries for the source term
+                for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+                    A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+
+                // restore the original state of the scv's volume variables
+                curVolVars = origVolVars;
+
+                // restore the current element solution
+                elemSol[0][pvIdx] = origPriVars[pvIdx];
+
+                // restore the undeflected state of the coupling context
+                this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalI, elemSol[0], pvIdx);
+            }
+
+            // restore original state of the flux vars cache in case of global caching.
+            // This has to be done in order to guarantee that everything is in an undeflected
+            // state before the assembly of another element is called. In the case of local caching
+            // this is obsolete because the elemFluxVarsCache used here goes out of scope after this.
+            // We only have to do this for the last primary variable, for all others the flux var cache
+            // is updated with the correct element volume variables before residual evaluations
+            if (enableGridFluxVarsCache)
+                gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
+        }
+
+        // analytic diff for source
+        else
+            this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, curElemVolVars, scv);
+
 
         // add flux derivatives for each scvf
         for (const auto& scvf : scvfs(fvGeometry))
@@ -859,8 +935,14 @@ public:
             return residual;
         }
         else
+        {
+            // compute potential additional derivatives caused by the coupling
+            const auto origResidual = this->evalLocalResidual();
+            this->couplingManager().evalAdditionalDomainDerivatives(domainI, *this, origResidual, A, gridVariables);
+
             // return element residual
-            return this->evalLocalResidual()[0];
+            return origResidual[0];
+        }
     }
 
     /*!
@@ -872,6 +954,42 @@ public:
     template<std::size_t otherId, class JacobianBlock, class GridVariables>
     void assembleJacobianCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
                                   const LocalResidualValues& res, GridVariables& gridVariables)
+    {
+        ///////////////////////////////////////////////////////
+        // Calculate derivatives of all dofs in the element  //
+        // with respect to all dofs in the coupling stencil. //
+        ///////////////////////////////////////////////////////
+        static const bool useNumDiff = getParam<bool>("Assembly.UseNumericDifferentiationForCouplingDerivatives", false);
+        if (useNumDiff)
+            assembleJacobianCouplingNumeric(domainJ, A, res, gridVariables);
+        else
+            this->couplingManager().addCouplingDerivatives(A, domainI, *this, domainJ);
+
+        if constexpr (Problem::enableInternalDirichletConstraints())
+        {
+            const auto globalI = this->fvGeometry().gridGeometry().elementMapper().index(this->element());
+            const auto& stencil = this->couplingManager().couplingStencil(domainI, this->element(), domainJ);
+            for (const auto globalJ : stencil)
+            {
+                const auto& scv = this->fvGeometry().scv(globalI);
+                const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
+                if (internalDirichletConstraints.any())
+                    for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+                        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+                            if (internalDirichletConstraints[eqIdx])
+                                A[globalI][globalJ][eqIdx][pvIdx] = 0.0;
+            }
+        }
+    }
+
+    /*!
+     * \brief Computes the derivatives with respect to the given element and adds them
+     *        to the global matrix.
+     * \todo This is just a copy of the assembleJacobianCoupling of the numeric-diff local assembler
+     */
+    template<std::size_t otherId, class JacobianBlock, class GridVariables>
+    void assembleJacobianCouplingNumeric(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
+                                         const LocalResidualValues& res, GridVariables& gridVariables)
     {
         ////////////////////////////////////////////////////////////////////////////////////////////////////////
         // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
@@ -882,33 +1000,89 @@ public:
         const auto& fvGeometry = this->fvGeometry();
         const auto& gridGeometry = fvGeometry.gridGeometry();
         auto&& curElemVolVars = this->curElemVolVars();
-        // auto&& elemFluxVarsCache = this->elemFluxVarsCache();
+        auto&& elemFluxVarsCache = this->elemFluxVarsCache();
 
         // get stencil informations
         const auto globalI = gridGeometry.elementMapper().index(element);
         const auto& stencil = this->couplingManager().couplingStencil(domainI, element, domainJ);
+        const auto& curSolJ = this->curSol()[domainJ];
 
-        for (const auto globalJ : stencil)
-        {
-            const auto& elementJ = this->assembler().gridGeometry(domainJ).element(globalJ);
-            this->couplingManager().addCouplingDerivatives(A[globalI][globalJ], domainI, element, fvGeometry, curElemVolVars, domainJ, elementJ);
+        // make sure the flux variables cache does not contain any artifacts from prior differentiation
+        elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
 
-            // add the current partial derivatives to the global jacobian matrix
-            if constexpr (Problem::enableInternalDirichletConstraints())
+        // convenience lambda for call to update self
+        auto updateCoupledVariables = [&] ()
+        {
+            // Update ourself after the context has been modified. Depending on the
+            // type of caching, other objects might have to be updated. All ifs can be optimized away.
+            if (enableGridFluxVarsCache)
             {
-                const auto& scv = fvGeometry.scv(globalI);
-                const auto internalDirichletConstraints = this->problem().hasInternalDirichletConstraint(this->element(), scv);
-                if (internalDirichletConstraints.any())
+                if (enableGridVolVarsCache)
                 {
-                    for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
-                        for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
-                            if (internalDirichletConstraints[eqIdx])
-                                A[globalI][globalJ][eqIdx][pvIdx] = 0.0;
+                    this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), gridVariables.gridFluxVarsCache());
+                    this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache);
                 }
+                else
+                {
+                    this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, gridVariables.gridFluxVarsCache());
+                    this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache);
+                }
+            }
+            else
+            {
+                if (enableGridVolVarsCache)
+                    this->couplingManager().updateCoupledVariables(domainI, *this, gridVariables.curGridVolVars(), elemFluxVarsCache);
+                else
+                    this->couplingManager().updateCoupledVariables(domainI, *this, curElemVolVars, elemFluxVarsCache);
+            }
+        };
+
+        for (const auto& globalJ : stencil)
+        {
+            // undeflected privars and privars to be deflected
+            const auto origPriVarsJ = curSolJ[globalJ];
+            auto priVarsJ = origPriVarsJ;
+
+            const auto origResidual = this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ)[0];
+
+            for (int pvIdx = 0; pvIdx < JacobianBlock::block_type::cols; ++pvIdx)
+            {
+                auto evalCouplingResidual = [&](Scalar priVar)
+                {
+                    // update the volume variables and the flux var cache
+                    priVarsJ[pvIdx] = priVar;
+                    this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
+                    updateCoupledVariables();
+                    return this->couplingManager().evalCouplingResidual(domainI, *this, domainJ, globalJ)[0];
+                };
+
+                // derive the residuals numerically
+                LocalResidualValues partialDeriv(0.0);
+                const auto& paramGroup = this->assembler().problem(domainJ).paramGroup();
+                static const int numDiffMethod = getParamFromGroup<int>(paramGroup, "Assembly.NumericDifferenceMethod");
+                static const auto epsCoupl = this->couplingManager().numericEpsilon(domainJ, paramGroup);
+                NumericDifferentiation::partialDerivative(evalCouplingResidual, priVarsJ[pvIdx], partialDeriv, origResidual,
+                                                          epsCoupl(priVarsJ[pvIdx], pvIdx), numDiffMethod);
+
+                // add the current partial derivatives to the global jacobian matrix
+                for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+                    A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
+
+                // restore the current element solution
+                priVarsJ[pvIdx] = origPriVarsJ[pvIdx];
+
+                // restore the undeflected state of the coupling context
+                this->couplingManager().updateCouplingContext(domainI, *this, domainJ, globalJ, priVarsJ, pvIdx);
             }
+
+            // restore original state of the flux vars cache and/or vol vars in case of global caching.
+            // This has to be done in order to guarantee that the undeflected residual computation done
+            // above is correct when jumping to the next coupled element.
+            updateCoupledVariables();
         }
     }
-};
+
+}; // implicit assebler with analytic Jacobian
 
 } // end namespace Dumux
 
diff --git a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
index 184fffe06813e17520a0dcb73076d3d5785f5255..7e016d799c5c6bc191551bcc32a05cabacd6f0b5 100644
--- a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
+++ b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
@@ -76,10 +76,10 @@ public:
                               const Problem& problem,
                               const Element& element,
                               const FVElementGeometry& fvGeometry,
-                              const VolumeVariables& curVolVars,
+                              const ElementVolumeVariables& curElemVolVars,
                               const SubControlVolume& scv) const
     {
-        problem.addSourceDerivatives(partialDerivatives, element, fvGeometry, curVolVars, scv);
+        problem.addSourceDerivatives(partialDerivatives, element, fvGeometry, curElemVolVars, scv);
     }
 
     //! Flux derivatives for the cell-centered tpfa scheme
diff --git a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
index 9475f5c797dd4e38e037364b9e9ae660668a986c..f8658eb5cd879c83493f9fa98f9e9bb563fdac40 100644
--- a/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
+++ b/dumux/porousmediumflow/2p/incompressiblelocalresidual.hh
@@ -128,7 +128,7 @@ public:
                               const Problem& problem,
                               const Element& element,
                               const FVElementGeometry& fvGeometry,
-                              const VolumeVariables& curVolVars,
+                              const ElementVolumeVariables& curElemVolVars,
                               const SubControlVolume& scv) const
     { /* TODO maybe forward to problem for the user to implement the source derivatives?*/ }
 
diff --git a/dumux/porousmediumflow/richards/localresidual.hh b/dumux/porousmediumflow/richards/localresidual.hh
index 5434833c98a58d31d9c8d3ec9c4e5ec0828ec194..31a121ae2d9b9ce1600d0a4720a2f15825691b15 100644
--- a/dumux/porousmediumflow/richards/localresidual.hh
+++ b/dumux/porousmediumflow/richards/localresidual.hh
@@ -224,9 +224,11 @@ public:
                               const Problem& problem,
                               const Element& element,
                               const FVElementGeometry& fvGeometry,
-                              const VolumeVariables& curVolVars,
+                              const ElementVolumeVariables& curElemVolVars,
                               const SubControlVolume& scv) const
-    { /* TODO maybe forward to problem for the user to implement the source derivatives?*/ }
+    {
+        problem.addSourceDerivatives(partialDerivatives, element, fvGeometry, curElemVolVars, scv);
+    }
 
     /*!
      * \brief Adds flux derivatives for wetting and non-wetting phase for cell-centered FVM using TPFA
diff --git a/dumux/porousmediumflow/tracer/localresidual.hh b/dumux/porousmediumflow/tracer/localresidual.hh
index 702ff39e990eaca7d947cbba46f90957aa684811..047ff12eff0a48bc504a12bd72aa56dd3b4a2501 100644
--- a/dumux/porousmediumflow/tracer/localresidual.hh
+++ b/dumux/porousmediumflow/tracer/localresidual.hh
@@ -229,7 +229,7 @@ public:
                               const Problem& problem,
                               const Element& element,
                               const FVElementGeometry& fvGeometry,
-                              const VolumeVariables& curVolVars,
+                              const ElementVolumeVariables& curElemVolVars,
                               const SubControlVolume& scv) const
     {
         // TODO maybe forward to the problem? -> necessary for reaction terms
diff --git a/test/multidomain/embedded/1d3d/1p_richards/CMakeLists.txt b/test/multidomain/embedded/1d3d/1p_richards/CMakeLists.txt
index 690280bb6651b77c5047f8eac78781ad31bbda0a..fadd9441b3d967acd598675a40457b42d0ec524d 100644
--- a/test/multidomain/embedded/1d3d/1p_richards/CMakeLists.txt
+++ b/test/multidomain/embedded/1d3d/1p_richards/CMakeLists.txt
@@ -2,6 +2,7 @@ dumux_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfatpfa
               LABELS multidomain multidomain_embedded 1p richards
               SOURCES main.cc
               COMPILE_DEFINITIONS SOILTYPETAG=SoilCC
+              COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::numeric
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMAKE_GUARD dune-foamgrid_FOUND
               TIMEOUT 1500
@@ -13,10 +14,27 @@ dumux_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfatpfa
                         --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfatpfa params.input \
                                    -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfatpfa")
 
-dumux_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfabox
+dune_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff
+              LABELS multidomain multidomain_embedded 1p richards
+              SOURCES main.cc
+              COMPILE_DEFINITIONS SOILTYPETAG=SoilCC
+              COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::analytic
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMAKE_GUARD dune-foamgrid_FOUND
+              TIMEOUT 1500
+              CMD_ARGS  --script fuzzy
+                        --files ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_1d-reference.vtp
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff_1d-00004.vtp
+                                ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_3d-reference.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff_3d-00004.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff params.input \
+                      -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfatpfa_anadiff")
+
+dune_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfabox
               LABELS multidomain multidomain_embedded 1p richards
               SOURCES main.cc
               COMPILE_DEFINITIONS SOILTYPETAG=SoilBox
+              COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::numeric
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMAKE_GUARD dune-foamgrid_FOUND
               TIMEOUT 1500
@@ -26,5 +44,55 @@ dumux_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfabox
                                 ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfabox_3d-reference.vtu
                                 ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox_3d-00004.vtu
                         --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox params.input \
-                                   -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfabox")
+                      -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfabox")
+
+dune_add_test(NAME test_md_embedded_1d3d_1p_richards_tpfabox_anadiff
+              LABELS multidomain multidomain_embedded 1p richards
+              SOURCES main.cc
+              COMPILE_DEFINITIONS SOILTYPETAG=SoilBox
+              COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::analytic
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMAKE_GUARD dune-foamgrid_FOUND
+              TIMEOUT 1500
+              CMD_ARGS  --script fuzzy
+                        --files ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_1d-reference.vtp
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox_anadiff_1d-00004.vtp
+                                ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfabox_3d-reference.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox_anadiff_3d-00004.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_tpfabox_anadiff params.input \
+                      -Vtk.OutputName test_md_embedded_1d3d_1p_richards_tpfabox_anadiff")
+
+dune_add_test(NAME test_md_embedded_1d3d_1p_richards_boxbox
+              LABELS multidomain multidomain_embedded 1p richards
+              SOURCES main.cc
+              COMPILE_DEFINITIONS SOILTYPETAG=SoilBox
+              COMPILE_DEFINITIONS ROOTTYPETAG=RootBox
+              COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::numeric
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMAKE_GUARD dune-foamgrid_FOUND
+              TIMEOUT 1500
+              CMD_ARGS  --script fuzzy
+                        --files ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_1d-reference.vtp
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_1d-00004.vtp
+                                ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_boxbox_3d-reference.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_3d-00004.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox params.input \
+                      -Vtk.OutputName test_md_embedded_1d3d_1p_richards_boxbox")
+
+dune_add_test(NAME test_md_embedded_1d3d_1p_richards_boxbox_anadiff
+              LABELS multidomain multidomain_embedded 1p richards
+              SOURCES main.cc
+              COMPILE_DEFINITIONS SOILTYPETAG=SoilBox
+              COMPILE_DEFINITIONS ROOTTYPETAG=RootBox
+              COMPILE_DEFINITIONS DIFFMETHOD=DiffMethod::analytic
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMAKE_GUARD dune-foamgrid_FOUND
+              TIMEOUT 1500
+              CMD_ARGS  --script fuzzy
+                        --files ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_1d-reference.vtp
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_anadiff_1d-00004.vtp
+                                ${CMAKE_SOURCE_DIR}/test/references/test_md_embedded_1d3d_1p_richards_boxbox_3d-reference.vtu
+                                ${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_anadiff_3d-00004.vtu
+                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_md_embedded_1d3d_1p_richards_boxbox_anadiff params.input \
+                      -Vtk.OutputName test_md_embedded_1d3d_1p_richards_boxbox_anadiff")
 dune_symlink_to_source_files(FILES "params.input")
diff --git a/test/multidomain/embedded/1d3d/1p_richards/main.cc b/test/multidomain/embedded/1d3d/1p_richards/main.cc
index dce774e34115370ad78830630388d507d2017924..4a99049f636f3d9c9cb4fb17202ec73973cd73e3 100644
--- a/test/multidomain/embedded/1d3d/1p_richards/main.cc
+++ b/test/multidomain/embedded/1d3d/1p_richards/main.cc
@@ -25,6 +25,13 @@
 
 #include <config.h>
 
+#ifndef SOILTYPETAG
+#define SOILTYPETAG SoilCC
+#endif
+#ifndef ROOTTYPETAG
+#define ROOTTYPETAG RootCC
+#endif
+
 #include <ctime>
 #include <iostream>
 #include <fstream>
@@ -56,12 +63,12 @@ namespace Properties {
 template<class TypeTag>
 struct CouplingManager<TypeTag, TTag::SOILTYPETAG>
 {
-    using Traits = MultiDomainTraits<TypeTag, Properties::TTag::Root>;
+    using Traits = MultiDomainTraits<TypeTag, Properties::TTag::ROOTTYPETAG>;
     using type = EmbeddedCouplingManager1d3d<Traits, EmbeddedCouplingMode::average>;
 };
 
 template<class TypeTag>
-struct CouplingManager<TypeTag, TTag::Root>
+struct CouplingManager<TypeTag, TTag::ROOTTYPETAG>
 {
     using Traits = MultiDomainTraits<Properties::TTag::SOILTYPETAG, TypeTag>;
     using type = EmbeddedCouplingManager1d3d<Traits, EmbeddedCouplingMode::average>;
@@ -70,11 +77,11 @@ struct CouplingManager<TypeTag, TTag::Root>
 template<class TypeTag>
 struct PointSource<TypeTag, TTag::SOILTYPETAG> { using type = typename GetPropType<TypeTag, Properties::CouplingManager>::PointSourceTraits::template PointSource<0>; };
 template<class TypeTag>
-struct PointSource<TypeTag, TTag::Root> { using type = typename GetPropType<TypeTag, Properties::CouplingManager>::PointSourceTraits::template PointSource<1>; };
+struct PointSource<TypeTag, TTag::ROOTTYPETAG> { using type = typename GetPropType<TypeTag, Properties::CouplingManager>::PointSourceTraits::template PointSource<1>; };
 template<class TypeTag>
 struct PointSourceHelper<TypeTag, TTag::SOILTYPETAG> { using type = typename GetPropType<TypeTag, Properties::CouplingManager>::PointSourceTraits::template PointSourceHelper<0>; };
 template<class TypeTag>
-struct PointSourceHelper<TypeTag, TTag::Root> { using type = typename GetPropType<TypeTag, Properties::CouplingManager>::PointSourceTraits::template PointSourceHelper<1>; };
+struct PointSourceHelper<TypeTag, TTag::ROOTTYPETAG> { using type = typename GetPropType<TypeTag, Properties::CouplingManager>::PointSourceTraits::template PointSourceHelper<1>; };
 
 } // end namespace Properties
 } // end namespace Dumux
@@ -95,7 +102,7 @@ int main(int argc, char** argv)
 
     // Define the sub problem type tags
     using BulkTypeTag = Properties::TTag::SOILTYPETAG;
-    using LowDimTypeTag = Properties::TTag::Root;
+    using LowDimTypeTag = Properties::TTag::ROOTTYPETAG;
 
     // try to create a grid (from the given grid file or the input file)
     // for both sub-domains
@@ -188,7 +195,7 @@ int main(int argc, char** argv)
     timeLoop->setMaxTimeStepSize(maxDt);
 
     // the assembler with time loop for instationary problem
-    using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+    using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DIFFMETHOD>;
     auto assembler = std::make_shared<Assembler>(std::make_tuple(bulkProblem, lowDimProblem),
                                                  std::make_tuple(bulkFvGridGeometry, lowDimFvGridGeometry),
                                                  std::make_tuple(bulkGridVariables, lowDimGridVariables),
diff --git a/test/multidomain/embedded/1d3d/1p_richards/params.input b/test/multidomain/embedded/1d3d/1p_richards/params.input
index fe09ff1a666d46b5cec8437e7599da4534f03b63..408a9000750c827380c0aa33cb2ba239363c8dee 100644
--- a/test/multidomain/embedded/1d3d/1p_richards/params.input
+++ b/test/multidomain/embedded/1d3d/1p_richards/params.input
@@ -35,3 +35,7 @@ Problem.Name = 1d
 InitialSoilPressure =  -0.9429e4 # [Pa]
 InitialRootPressure = -1.2e6 # [Pa]
 TranspirationRate = 2.15e-8 # [kg / s]
+
+[Assembly]
+UseNumericDifferentiationForSourceDerivatives = true
+UseNumericDifferentiationForCouplingDerivatives = true
diff --git a/test/multidomain/embedded/1d3d/1p_richards/problem_root.hh b/test/multidomain/embedded/1d3d/1p_richards/problem_root.hh
index fd2cf0592801e04e4316ccf8106124cf0b24ecb9..801137810c561e2b27d148ceec8cf23e4bdcaddc 100644
--- a/test/multidomain/embedded/1d3d/1p_richards/problem_root.hh
+++ b/test/multidomain/embedded/1d3d/1p_richards/problem_root.hh
@@ -32,6 +32,7 @@
 #include <dumux/common/parameters.hh>
 #include <dumux/common/properties.hh>
 #include <dumux/discretization/cctpfa.hh>
+#include <dumux/discretization/box.hh>
 
 #include <dumux/porousmediumflow/1p/model.hh>
 #include <dumux/porousmediumflow/problem.hh>
@@ -50,7 +51,9 @@ namespace Properties {
 
 // Create new type tags
 namespace TTag {
-struct Root { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
+struct Root { using InheritsFrom = std::tuple<OneP>; };
+struct RootCC { using InheritsFrom = std::tuple<Root, CCTpfaModel>; };
+struct RootBox { using InheritsFrom = std::tuple<Root, BoxModel>; };
 } // end namespace TTag
 
 // Set the grid type