diff --git a/ebos/eclbasevanguard.hh b/ebos/eclbasevanguard.hh
index b6d2f98b8..e06b580fd 100644
--- a/ebos/eclbasevanguard.hh
+++ b/ebos/eclbasevanguard.hh
@@ -800,7 +800,7 @@ protected:
             const unsigned int elemIdx = elemMapper.index(element);
             cellCenterDepth_[elemIdx] = cellCenterDepth(element);
 
-            if (this->eclState_->aquifer().hasNumericalAquifer()) {
+            if (!num_aqu_cells.empty()) {
                 const unsigned int global_index = cartesianIndex(elemIdx);
                 const auto search = num_aqu_cells.find(global_index);
                 if (search != num_aqu_cells.end()) {
diff --git a/ebos/eclcpgridvanguard.hh b/ebos/eclcpgridvanguard.hh
index 78a826542..8a8e296b1 100644
--- a/ebos/eclcpgridvanguard.hh
+++ b/ebos/eclcpgridvanguard.hh
@@ -341,7 +341,6 @@ protected:
     void createGrids_()
     {
         grid_.reset(new Dune::CpGrid());
-
         grid_->processEclipseFormat(mpiRank == 0 ? &this->eclState().getInputGrid()
                                                  : nullptr,
                                     /*isPeriodic=*/false,
diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh
index c24d78fe8..e7d03b6ef 100644
--- a/ebos/equil/initstateequil.hh
+++ b/ebos/equil/initstateequil.hh
@@ -1935,10 +1935,8 @@ private:
         using CellPos = typename PhaseSat::Position;
         using CellID  = std::remove_cv_t<std::remove_reference_t<
             decltype(std::declval<CellPos>().cell)>>;
-        // TODO: might not needed
-        ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
         const auto num_aqu_cells = aquifer.allAquiferCells();
-        this->cellLoop(cells, [this, &eqreg,  &ptable, &elemMapper, &num_aqu_cells, &psat]
+        this->cellLoop(cells, [this, &eqreg,  &ptable, &num_aqu_cells, &psat]
             (const CellID                 cell,
              Details::PhaseQuantityValue& pressures,
              Details::PhaseQuantityValue& saturations,
@@ -1964,7 +1962,6 @@ private:
                     OpmLog::info(msg);
 
                     if (aqu_cell->init_pressure) {
-                        // TODO: NOT totally sure what we should do here to employ the pressure specified by AQUNUM
                         const double pres = *(aqu_cell->init_pressure);
                         pressures = {pres, pres, pres};
                     }
diff --git a/opm/simulators/aquifers/AquiferNumerical.hpp b/opm/simulators/aquifers/AquiferNumerical.hpp
index fd34c0c2d..c88b08879 100644
--- a/opm/simulators/aquifers/AquiferNumerical.hpp
+++ b/opm/simulators/aquifers/AquiferNumerical.hpp
@@ -22,8 +22,7 @@
 #define OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
 
 #include <opm/output/data/Aquifer.hpp>
-
-#include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.hpp>
+#include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
 
 namespace Opm
 {
@@ -53,7 +52,7 @@ public:
                      const std::unordered_map<int, int>& cartesian_to_compressed,
                      const Simulator& ebos_simulator,
                      const int* global_cell)
-    : aquifer_(aquifer)
+    : id_(aquifer.id())
     , ebos_simulator_(ebos_simulator)
     , flux_rate_(0.)
     , cumulative_flux_(0.)
@@ -85,7 +84,7 @@ public:
     Opm::data::AquiferData aquiferData() const
     {
         data::AquiferData data;
-        data.aquiferID = this->aquifer_.id();
+        data.aquiferID = this->id_;
         data.initPressure = this->init_pressure_;
         data.pressure = this->pressure_;
         data.fluxRate = this->flux_rate_;
@@ -103,7 +102,7 @@ public:
     }
 
 private:
-    const Opm::SingleNumericalAquifer& aquifer_;
+    const size_t id_;
     const Simulator& ebos_simulator_;
     double flux_rate_; // aquifer influx rate
     double cumulative_flux_; // cumulative aquifer influx
@@ -144,7 +143,6 @@ private:
             sum_pv += pv;
         }
 
-        assert(sum_pv > 0.);
         return sum_pv_pressure/ sum_pv;
     }
 
diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
index 7be5424bc..7fb8a3d4a 100644
--- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
+++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
@@ -71,6 +71,11 @@ BlackoilAquiferModel<TypeTag>::initFromRestart(const std::vector<data::AquiferDa
             aquifer.initFromRestart(aquiferSoln);
         }
     }
+    if (aquiferNumericalActive()) {
+        for (auto& aquifer : this->aquifers_numerical) {
+            aquifer.initFromRestart(aquiferSoln);
+        }
+    }
 }
 
 template <typename TypeTag>
@@ -195,14 +200,13 @@ BlackoilAquiferModel<TypeTag>::init()
     }
 
     if (aquifer.hasNumericalAquifer()) {
-        const auto aquifers = aquifer.numericalAquifers().aquifers();
+        const auto& num_aquifers = aquifer.numericalAquifers().aquifers();
         const auto& ugrid = simulator_.vanguard().grid();
         const int number_of_cells = simulator_.gridView().size(0);
         const int* global_cell = Opm::UgGridHelpers::globalCell(ugrid);
         const std::unordered_map<int, int> cartesian_to_compressed = cartesianToCompressed(number_of_cells,
                                                                                            global_cell);
-        // for (const auto& elem : ers().aquifers()) {
-        for (const auto& [id, aqu] : aquifers) {
+        for ([[maybe_unused]]const auto& [id, aqu] : num_aquifers) {
             this->aquifers_numerical.emplace_back(aqu,
                   cartesian_to_compressed, this->simulator_, global_cell);
         }
@@ -210,12 +214,6 @@ BlackoilAquiferModel<TypeTag>::init()
 }
 template <typename TypeTag>
 bool
-BlackoilAquiferModel<TypeTag>::aquiferActive() const
-{
-    return (aquiferCarterTracyActive() || aquiferFetkovichActive());
-}
-template <typename TypeTag>
-bool
 BlackoilAquiferModel<TypeTag>::aquiferCarterTracyActive() const
 {
     return !aquifers_CarterTracy.empty();