diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index 500c31102..418be36a6 100644
--- a/CMakeLists_files.cmake
+++ b/CMakeLists_files.cmake
@@ -161,6 +161,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/simulators/aquifers/AquiferInterface.hpp
   opm/simulators/aquifers/AquiferCarterTracy.hpp
   opm/simulators/aquifers/AquiferFetkovich.hpp
+  opm/simulators/aquifers/AquiferNumerical.hpp
   opm/simulators/aquifers/BlackoilAquiferModel.hpp
   opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
   opm/simulators/linalg/bda/BdaBridge.hpp
diff --git a/compareECLFiles.cmake b/compareECLFiles.cmake
index 296141376..e60ce003c 100755
--- a/compareECLFiles.cmake
+++ b/compareECLFiles.cmake
@@ -283,6 +283,22 @@ add_test_compareECLFiles(CASENAME fetkovich_2d
                          REL_TOL ${rel_tol}
                          DIR aquifer-fetkovich)
 
+add_test_compareECLFiles(CASENAME numerical_aquifer_3d_2aqu
+                         FILENAME 3D_2AQU_NUM
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${rel_tol}
+                         DIR aquifer-num
+                         TEST_ARGS --relaxed-max-pv-fraction=0 --tolerance-cnv=0.00003 --time-step-control=pid --linsolver=cpr)
+
+add_test_compareECLFiles(CASENAME numerical_aquifer_3d_1aqu
+                         FILENAME 3D_1AQU_3CELLS
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${rel_tol}
+                         DIR aquifer-num
+                         TEST_ARGS --relaxed-max-pv-fraction=0 --tolerance-cnv=0.00003 --time-step-control=pid --linsolver=cpr)
+
 add_test_compareECLFiles(CASENAME spe3
                          FILENAME SPE3CASE1
                          SIMULATOR flow
@@ -998,4 +1014,20 @@ if(MPI_FOUND)
                                        REL_TOL ${rel_tol_parallel}
                                        DIR parallel_fieldprops
                                        TEST_ARGS --linear-solver-reduction=1e-7 --tolerance-cnv=5e-6 --tolerance-mb=1e-6)
+
+  add_test_compare_parallel_simulation(CASENAME numerical_aquifer_3d_2aqu
+                                       FILENAME 3D_2AQU_NUM
+                                       SIMULATOR flow
+                                       ABS_TOL 0.12
+                                       REL_TOL ${coarse_rel_tol_parallel}
+                                       DIR aquifer-num
+                                       TEST_ARGS --relaxed-max-pv-fraction=0 --tolerance-cnv=0.00003 --time-step-control=pid --linsolver=cpr)
+
+  add_test_compare_parallel_simulation(CASENAME numerical_aquifer_3d_1aqu
+                                       FILENAME 3D_1AQU_3CELLS
+                                       SIMULATOR flow
+                                       ABS_TOL ${abs_tol_parallel}
+                                       REL_TOL ${coarse_rel_tol_parallel}
+                                       DIR aquifer-num
+                                       TEST_ARGS --relaxed-max-pv-fraction=0 --tolerance-cnv=0.00003 --time-step-control=pid --linsolver=cpr)
 endif()
diff --git a/ebos/eclbasevanguard.hh b/ebos/eclbasevanguard.hh
index 7035f9d6a..2d3f2cada 100644
--- a/ebos/eclbasevanguard.hh
+++ b/ebos/eclbasevanguard.hh
@@ -210,6 +210,7 @@ public:
 protected:
     static const int dimension = Grid::dimension;
     using Element = typename GridView::template Codim<0>::Entity;
+    using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
 
 
 public:
@@ -637,6 +638,9 @@ public:
     const std::string& caseName() const
     { return caseName_; }
 
+    const CartesianIndexMapper& cartesianMapper() const
+    {  return asImp_().cartesianIndexMapper(); }
+
     /*!
      * \brief Returns the number of logically Cartesian cells in each direction
      */
@@ -787,10 +791,22 @@ protected:
         ElementMapper elemMapper(this->gridView(), Dune::mcmgElementLayout());
         auto elemIt = this->gridView().template begin</*codim=*/0>();
         const auto& elemEndIt = this->gridView().template end</*codim=*/0>();
+
+        const auto num_aqu_cells = this->eclState_->aquifer().numericalAquifers().allAquiferCells();
+
         for (; elemIt != elemEndIt; ++elemIt) {
             const Element& element = *elemIt;
             const unsigned int elemIdx = elemMapper.index(element);
             cellCenterDepth_[elemIdx] = cellCenterDepth(element);
+
+            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()) {
+                    // updating the cell depth using aquifer cell depth
+                    cellCenterDepth_[elemIdx] = search->second->depth;
+                }
+            }
         }
     }
 
diff --git a/ebos/eclequilinitializer.hh b/ebos/eclequilinitializer.hh
index 5c940f5c2..a5aef79cf 100644
--- a/ebos/eclequilinitializer.hh
+++ b/ebos/eclequilinitializer.hh
@@ -101,6 +101,7 @@ public:
         EQUIL::DeckDependent::InitialStateComputer<TypeTag> initialState(materialLawManager,
                                                                          eclState,
                                                                          vanguard.gridView(),
+                                                                         vanguard.cartesianMapper(),
                                                                          simulator.problem().gravity()[dimWorld - 1]);
 
         // copy the result into the array of initial fluid states
diff --git a/ebos/ecltransmissibility.hh b/ebos/ecltransmissibility.hh
index 9d2bcaf3b..1d77cb486 100644
--- a/ebos/ecltransmissibility.hh
+++ b/ebos/ecltransmissibility.hh
@@ -805,7 +805,7 @@ private:
                 // Discard the NNC if it is between active cell and inactive cell
                 std::ostringstream sstr;
                 sstr << "NNC between active and inactive cells ("
-                     << low << " -> " << high << ")";
+                     << low << " -> " << high << ") with globalcell is (" << c1 << "->" << c2 <<")";
                 Opm::OpmLog::warning(sstr.str());
                 continue;
             }
diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh
index c5eca90cb..c867a8742 100644
--- a/ebos/equil/initstateequil.hh
+++ b/ebos/equil/initstateequil.hh
@@ -48,6 +48,7 @@
 #include <opm/parser/eclipse/EclipseState/Tables/PdvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SaltvdTable.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
+#include <fmt/format.h>
 
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 #include <opm/material/fluidstates/SimpleModularFluidState.hpp>
@@ -62,6 +63,7 @@
 #include <type_traits>
 #include <utility>
 #include <vector>
+#include <string>
 
 namespace Opm {
 
@@ -1566,12 +1568,14 @@ class InitialStateComputer
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
     using Grid = GetPropType<TypeTag, Properties::Grid>;
+    using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
 
 public:
     template<class MaterialLawManager>
     InitialStateComputer(MaterialLawManager& materialLawManager,
                          const Opm::EclipseState& eclipseState,
                          const GridView& gridView,
+                         const CartesianIndexMapper& cartMapper,
                          const double grav = Opm::unit::gravity,
                          const bool applySwatInit = true)
         : temperature_(gridView.size(/*codim=*/0)),
@@ -1581,7 +1585,8 @@ public:
           sat_(FluidSystem::numPhases,
           std::vector<double>(gridView.size(/*codim=*/0))),
           rs_(gridView.size(/*codim=*/0)),
-          rv_(gridView.size(/*codim=*/0))
+          rv_(gridView.size(/*codim=*/0)),
+          cartesianIndexMapper_(cartMapper)
     {
         //Check for presence of kw SWATINIT
         if (applySwatInit) {
@@ -1590,8 +1595,10 @@ public:
             }
         }
 
-        //Querry cell depth, cell top-bottom.
-        updateCellProps_(gridView);
+        // Querry cell depth, cell top-bottom.
+        // numerical aquifer cells might be specified with different depths.
+        const auto& num_aquifers = eclipseState.aquifer().numericalAquifers();
+        updateCellProps_(gridView, num_aquifers);
 
         // Get the equilibration records.
         const std::vector<Opm::EquilRecord> rec = getEquil(eclipseState);
@@ -1708,6 +1715,9 @@ public:
         const auto& comm = gridView.comm();
         calcPressSatRsRv(eqlmap, rec, materialLawManager, comm, grav);
 
+        // modify the pressure and saturation for numerical aquifer cells
+        applyNumericalAquifers_(gridView, num_aquifers);
+
         // Modify oil pressure in no-oil regions so that the pressures of present phases can
         // be recovered from the oil pressure and capillary relations.
     }
@@ -1769,12 +1779,14 @@ private:
     PVec sat_;
     Vec rs_;
     Vec rv_;
+    const CartesianIndexMapper& cartesianIndexMapper_;
     Vec swatInit_;
     Vec cellCenterDepth_;
     std::vector<std::pair<double,double>> cellZSpan_;
     std::vector<std::pair<double,double>> cellZMinMax_;
 
-    void updateCellProps_(const GridView& gridView)
+    void updateCellProps_(const GridView& gridView,
+                          const NumericalAquifers& aquifer)
     {
         ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
         int numElements = gridView.size(/*codim=*/0);
@@ -1784,15 +1796,77 @@ private:
 
         auto elemIt = gridView.template begin</*codim=*/0>();
         const auto& elemEndIt = gridView.template end</*codim=*/0>();
+        const auto num_aqu_cells = aquifer.allAquiferCells();
         for (; elemIt != elemEndIt; ++elemIt) {
             const Element& element = *elemIt;
             const unsigned int elemIdx = elemMapper.index(element);
             cellCenterDepth_[elemIdx] = Details::cellCenterDepth(element);
+            const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
+            if (!num_aqu_cells.empty()) {
+                const auto search = num_aqu_cells.find(cartIx);
+                if (search != num_aqu_cells.end()) {
+                    const auto* aqu_cell = num_aqu_cells.at(cartIx);
+                    cellCenterDepth_[elemIdx] = aqu_cell->depth;
+                }
+            }
             cellZSpan_[elemIdx] = Details::cellZSpan(element);
             cellZMinMax_[elemIdx] = Details::cellZMinMax(element);
         }
     }
 
+    void applyNumericalAquifers_(const GridView& gridView,
+                                 const NumericalAquifers& aquifer)
+    {
+        const auto num_aqu_cells = aquifer.allAquiferCells();
+        if (num_aqu_cells.empty()) return;
+
+        ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
+        auto elemIt = gridView.template begin</*codim=*/0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const Element& element = *elemIt;
+            const unsigned int elemIdx = elemMapper.index(element);
+            const auto cartIx = cartesianIndexMapper_.cartesianIndex(elemIdx);
+            const auto search = num_aqu_cells.find(cartIx);
+            if (search != num_aqu_cells.end()) {
+                // numerical aquifer cells are filled with water initially
+                const auto watPos = FluidSystem::waterPhaseIdx;
+                if (FluidSystem::phaseIsActive(watPos)) {
+                    this->sat_[watPos][elemIdx] = 1.;
+                } else {
+                    throw std::logic_error  { "Water phase has to be active for numerical aquifer case" };
+                }
+
+                const auto oilPos = FluidSystem::oilPhaseIdx;
+                if (FluidSystem::phaseIsActive(oilPos)) {
+                    this->sat_[oilPos][elemIdx] = 0.;
+                }
+
+                const auto gasPos = FluidSystem::gasPhaseIdx;
+                if (FluidSystem::phaseIsActive(gasPos)) {
+                    this->sat_[gasPos][elemIdx] = 0.;
+                }
+                const auto* aqu_cell = num_aqu_cells.at(cartIx);
+                const auto msg = fmt::format("FOR AQUIFER CELL AT ({}, {}, {}) OF NUMERICAL "
+                                             "AQUIFER {}, WATER SATURATION IS SET TO BE UNITY",
+                                             aqu_cell->I+1, aqu_cell->J+1, aqu_cell->K+1, aqu_cell->aquifer_id);
+                OpmLog::info(msg);
+
+                // if pressure is specified for numerical aquifers, we use these pressure values
+                // for numerical aquifer cells
+                if (aqu_cell->init_pressure) {
+                    const double pres = *(aqu_cell->init_pressure);
+                    this->pp_[watPos][elemIdx] = pres;
+                    if (FluidSystem::phaseIsActive(gasPos)) {
+                        this->pp_[gasPos][elemIdx] = pres;
+                    }
+                    if (FluidSystem::phaseIsActive(oilPos)) {
+                        this->pp_[oilPos][elemIdx] = pres;
+                    }
+                }
+            }
+        }
+    }
     template<class RMap>
     void setRegionPvtIdx(const Opm::EclipseState& eclState, const RMap& reg)
     {
@@ -1806,7 +1880,7 @@ private:
 
     template <class RMap, class MaterialLawManager, class Comm>
     void calcPressSatRsRv(const RMap& reg,
-                          const std::vector< Opm::EquilRecord >& rec,
+                          const std::vector<Opm::EquilRecord>& rec,
                           MaterialLawManager& materialLawManager,
                           const Comm& comm,
                           const double grav)
@@ -1905,16 +1979,15 @@ private:
     }
 
     template <class CellRange, class PressTable, class PhaseSat>
-    void equilibrateCellCentres(const CellRange&  cells,
-                                const EquilReg&   eqreg,
-                                const PressTable& ptable,
-                                PhaseSat&         psat)
+    void equilibrateCellCentres(const CellRange&         cells,
+                                const EquilReg&          eqreg,
+                                const PressTable&        ptable,
+                                PhaseSat&                psat)
     {
         using CellPos = typename PhaseSat::Position;
         using CellID  = std::remove_cv_t<std::remove_reference_t<
             decltype(std::declval<CellPos>().cell)>>;
-
-        this->cellLoop(cells, [this, &eqreg, &ptable, &psat]
+        this->cellLoop(cells, [this, &eqreg,  &ptable, &psat]
             (const CellID                 cell,
              Details::PhaseQuantityValue& pressures,
              Details::PhaseQuantityValue& saturations,
diff --git a/opm/simulators/aquifers/AquiferNumerical.hpp b/opm/simulators/aquifers/AquiferNumerical.hpp
new file mode 100644
index 000000000..d671508df
--- /dev/null
+++ b/opm/simulators/aquifers/AquiferNumerical.hpp
@@ -0,0 +1,228 @@
+/*
+  Copyright (C) 2020 Equinor ASA
+  Copyright (C) 2020 SINTEF Digital
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM 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.
+
+  OPM 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 OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
+#define OPM_AQUIFERNUMERICAL_HEADER_INCLUDED
+
+#include <opm/output/data/Aquifer.hpp>
+#include <opm/parser/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
+
+namespace Opm
+{
+template <typename TypeTag>
+class AquiferNumerical
+{
+public:
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
+
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
+
+    enum { dimWorld = GridView::dimensionworld };
+
+    static const auto waterPhaseIdx = FluidSystem::waterPhaseIdx;
+
+    static const int numEq = BlackoilIndices::numEq;
+
+    using Eval =  DenseAd::Evaluation<double, numEq>;
+    using Toolbox = Opm::MathToolbox<Eval>;
+
+    // Constructor
+    AquiferNumerical(const SingleNumericalAquifer& aquifer,
+                     const std::unordered_map<int, int>& cartesian_to_compressed,
+                     const Simulator& ebos_simulator,
+                     const int* global_cell)
+    : id_(aquifer.id())
+    , ebos_simulator_(ebos_simulator)
+    , flux_rate_(0.)
+    , cumulative_flux_(0.)
+    , global_cell_(global_cell)
+    {
+        this->cell_to_aquifer_cell_idx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
+
+        for (size_t idx = 0; idx < aquifer.numCells(); ++idx) {
+            const auto& cell = *(aquifer.getCellPrt(idx));
+            const int global_idx = cell.global_index;
+            const auto search = cartesian_to_compressed.find(global_idx);
+            // Due to parallelisation, the cell might not exist in the current process
+            if (search != cartesian_to_compressed.end()) {
+                const int cell_idx = cartesian_to_compressed.at(global_idx);
+                this->cell_to_aquifer_cell_idx_[cell_idx] = idx;
+            }
+        }
+    }
+
+    void initFromRestart([[maybe_unused]]const std::vector<data::AquiferData>& aquiferSoln)
+    {
+        // NOT handling Restart for now
+    }
+
+    void endTimeStep()
+    {
+        this->pressure_ = this->calculateAquiferPressure();
+        this->flux_rate_ = this->calculateAquiferFluxRate();
+        this->cumulative_flux_ += this->flux_rate_ * this->ebos_simulator_.timeStepSize();
+    }
+
+    Opm::data::AquiferData aquiferData() const
+    {
+        data::AquiferData data;
+        data.aquiferID = this->id_;
+        data.initPressure = this->init_pressure_;
+        data.pressure = this->pressure_;
+        data.fluxRate = this->flux_rate_;
+        data.volume = this->cumulative_flux_;
+        data.type = Opm::data::AquiferType::Numerical;
+        return data;
+    }
+
+    void initialSolutionApplied()
+    {
+        this->init_pressure_ = this->calculateAquiferPressure();
+        this->pressure_ = this->init_pressure_;
+        this->flux_rate_ = 0.;
+        this->cumulative_flux_ = 0.;
+    }
+
+private:
+    const size_t id_;
+    const Simulator& ebos_simulator_;
+    double flux_rate_; // aquifer influx rate
+    double cumulative_flux_; // cumulative aquifer influx
+    const int* global_cell_; // mapping to global index
+    double init_pressure_;
+    double pressure_; // aquifer pressure
+
+    // TODO: maybe unordered_map can also do the work to save memory?
+    std::vector<int> cell_to_aquifer_cell_idx_;
+
+    double calculateAquiferPressure() const
+    {
+        double sum_pressure_watervolume = 0.;
+        double sum_watervolume = 0.;
+
+        ElementContext  elem_ctx(this->ebos_simulator_);
+        const auto& gridView = this->ebos_simulator_.gridView();
+        auto elemIt = gridView.template begin</*codim=*/0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const auto& elem = *elemIt;
+            if (elem.partitionType() != Dune::InteriorEntity) {
+                continue;
+            }
+            elem_ctx.updatePrimaryStencil(elem);
+
+            const size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+            const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
+            if (idx < 0) {
+                continue;
+            }
+
+            elem_ctx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+            const auto& iq0 = elem_ctx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
+            const auto& fs = iq0.fluidState();
+
+            // TODO: the porosity of the cells are still wrong for numerical aquifer cells
+            // Because the dofVolume still based on the grid information.
+            // The pore volume is correct. Extra efforts will be done to get sensible porosity value here later.
+            const double water_saturation = fs.saturation(waterPhaseIdx).value();
+            const double porosity = iq0.porosity().value();
+            const double volume = elem_ctx.dofTotalVolume(0, 0);
+            // TODO: not sure we should use water pressure here
+            const double water_pressure_reservoir = fs.pressure(waterPhaseIdx).value();
+            const double water_volume = volume * porosity * water_saturation;
+            sum_pressure_watervolume += water_volume * water_pressure_reservoir;
+            sum_watervolume += water_volume;
+        }
+
+        const auto& comm = this->ebos_simulator_.vanguard().grid().comm();
+        comm.sum(&sum_pressure_watervolume, 1);
+        comm.sum(&sum_watervolume, 1);
+        return sum_pressure_watervolume / sum_watervolume;
+    }
+
+    double calculateAquiferFluxRate() const
+    {
+        double aquifer_flux = 0.;
+
+        ElementContext  elem_ctx(this->ebos_simulator_);
+        const auto& gridView = this->ebos_simulator_.gridView();
+        auto elemIt = gridView.template begin</*codim=*/0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const auto &elem = *elemIt;
+            if (elem.partitionType() != Dune::InteriorEntity) {
+                continue;
+            }
+            // elem_ctx.updatePrimaryStencil(elem);
+            elem_ctx.updateStencil(elem);
+
+            const size_t cell_index = elem_ctx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+            const int idx = this->cell_to_aquifer_cell_idx_[cell_index];
+            // we only need the first aquifer cell
+            if (idx != 0) {
+                continue;
+            }
+            elem_ctx.updateAllIntensiveQuantities();
+            elem_ctx.updateAllExtensiveQuantities();
+
+            const size_t num_interior_faces = elem_ctx.numInteriorFaces(/*timeIdx*/ 0);
+            // const auto &problem = elem_ctx.problem();
+            const auto &stencil = elem_ctx.stencil(0);
+            // const auto& inQuants = elem_ctx.intensiveQuantities(0, /*timeIdx*/ 0);
+
+            for (size_t face_idx = 0; face_idx < num_interior_faces; ++face_idx) {
+                const auto &face = stencil.interiorFace(face_idx);
+                // dof index
+                const size_t i = face.interiorIndex();
+                const size_t j = face.exteriorIndex();
+                // compressed index
+                // const size_t I = stencil.globalSpaceIndex(i);
+                const size_t J = stencil.globalSpaceIndex(j);
+
+                assert(stencil.globalSpaceIndex(i) == cell_index);
+
+                // we do not consider the flux within aquifer cells
+                // we only need the flux to the connections
+                if (this->cell_to_aquifer_cell_idx_[J] > 0) {
+                    continue;
+                }
+                const auto &exQuants = elem_ctx.extensiveQuantities(face_idx, /*timeIdx*/ 0);
+                const double water_flux = Toolbox::value(exQuants.volumeFlux(waterPhaseIdx));
+
+                const size_t up_id = water_flux >= 0. ? i : j;
+                const auto &intQuantsIn = elem_ctx.intensiveQuantities(up_id, 0);
+                const double invB = Toolbox::value(intQuantsIn.fluidState().invB(waterPhaseIdx));
+                const double face_area = face.area();
+                aquifer_flux += water_flux * invB * face_area;
+            }
+
+            // we only need to handle the first aquifer cell, we can exit loop here
+            break;
+        }
+
+        return aquifer_flux;
+    }
+};
+} // namespace Opm
+#endif
diff --git a/opm/simulators/aquifers/BlackoilAquiferModel.hpp b/opm/simulators/aquifers/BlackoilAquiferModel.hpp
index 17a87f7c8..af80516e6 100644
--- a/opm/simulators/aquifers/BlackoilAquiferModel.hpp
+++ b/opm/simulators/aquifers/BlackoilAquiferModel.hpp
@@ -34,6 +34,7 @@
 
 #include <opm/simulators/aquifers/AquiferCarterTracy.hpp>
 #include <opm/simulators/aquifers/AquiferFetkovich.hpp>
+#include <opm/simulators/aquifers/AquiferNumerical.hpp>
 
 #include <opm/grid/CpGrid.hpp>
 #include <opm/grid/polyhedralgrid.hh>
@@ -120,6 +121,7 @@ protected:
     // they share the base class
     mutable std::vector<AquiferCarterTracy_object> aquifers_CarterTracy;
     mutable std::vector<AquiferFetkovich_object> aquifers_Fetkovich;
+    std::vector<AquiferNumerical<TypeTag>> aquifers_numerical;
 
     // This initialization function is used to connect the parser objects with the ones needed by AquiferCarterTracy
     void init();
@@ -127,6 +129,7 @@ protected:
     bool aquiferActive() const;
     bool aquiferCarterTracyActive() const;
     bool aquiferFetkovichActive() const;
+    bool aquiferNumericalActive() const;
 };
 
 
diff --git a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
index 847fe545e..dc16bb615 100644
--- a/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
+++ b/opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp
@@ -18,6 +18,9 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <opm/grid/utility/cartesianToCompressed.hpp>
+#include "BlackoilAquiferModel.hpp"
+
 namespace Opm
 {
 
@@ -46,6 +49,12 @@ BlackoilAquiferModel<TypeTag>::initialSolutionApplied()
             aquifer.initialSolutionApplied();
         }
     }
+
+    if (this->aquiferNumericalActive()) {
+        for (auto& aquifer : this->aquifers_numerical) {
+            aquifer.initialSolutionApplied();
+        }
+    }
 }
 
 template <typename TypeTag>
@@ -62,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>
@@ -132,6 +146,12 @@ BlackoilAquiferModel<TypeTag>::endTimeStep()
             aquifer.endTimeStep();
         }
     }
+    if (aquiferNumericalActive()) {
+        for (auto& aquifer : this->aquifers_numerical) {
+            aquifer.endTimeStep();
+            this->simulator_.vanguard().grid().comm().barrier();
+        }
+    }
 }
 template <typename TypeTag>
 void
@@ -168,7 +188,7 @@ BlackoilAquiferModel<TypeTag>::init()
         return;
     }
 
-    // Get all the carter tracy aquifer properties data and put it in aquifers vector
+
     const auto& connections = aquifer.connections();
     for (const auto& aq : aquifer.ct()) {
         aquifers_CarterTracy.emplace_back(connections[aq.aquiferID],
@@ -179,12 +199,19 @@ BlackoilAquiferModel<TypeTag>::init()
         aquifers_Fetkovich.emplace_back(connections[aq.aquiferID],
                                         this->simulator_, aq);
     }
-}
-template <typename TypeTag>
-bool
-BlackoilAquiferModel<TypeTag>::aquiferActive() const
-{
-    return (aquiferCarterTracyActive() || aquiferFetkovichActive());
+
+    if (aquifer.hasNumericalAquifer()) {
+        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 ([[maybe_unused]]const auto& [id, aqu] : num_aquifers) {
+            this->aquifers_numerical.emplace_back(aqu,
+                  cartesian_to_compressed, this->simulator_, global_cell);
+        }
+    }
 }
 template <typename TypeTag>
 bool
@@ -199,6 +226,13 @@ BlackoilAquiferModel<TypeTag>::aquiferFetkovichActive() const
     return !aquifers_Fetkovich.empty();
 }
 
+template<typename TypeTag>
+bool
+BlackoilAquiferModel<TypeTag>::aquiferNumericalActive() const
+{
+    return !(this->aquifers_numerical.empty());
+}
+
 template<typename TypeTag>
 Opm::data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const {
     Opm::data::Aquifers data;
@@ -215,6 +249,14 @@ Opm::data::Aquifers BlackoilAquiferModel<TypeTag>::aquiferData() const {
             data[aqu_data.aquiferID] = aqu_data;
         }
     }
+
+    if (this->aquiferNumericalActive()) {
+        for (const auto& aqu : this->aquifers_numerical) {
+            Opm::data::AquiferData aqu_data = aqu.aquiferData();
+            data[aqu_data.aquiferID] = aqu_data;
+        }
+    }
+
     return data;
 }
 } // namespace Opm
diff --git a/tests/test_equil.cc b/tests/test_equil.cc
index b2bd9752b..f7a6b66da 100644
--- a/tests/test_equil.cc
+++ b/tests/test_equil.cc
@@ -475,7 +475,9 @@ BOOST_AUTO_TEST_CASE(DeckAllDead)
     Opm::GridManager gm(eclipseState.getInputGrid());
     const UnstructuredGrid& grid = *(gm.c_grid());
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 10.0);
+    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
+                                                                  eclipseState, simulator->vanguard().gridView(),
+                                                                  simulator->vanguard().cartesianMapper(), 10.0);
     const auto& pressures = comp.press();
     BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
     BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@@ -552,7 +554,9 @@ BOOST_AUTO_TEST_CASE(DeckWithCapillary)
     Opm::GridManager gm(eclipseState.getInputGrid());
     const UnstructuredGrid& grid = *(gm.c_grid());
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 10.0);
+    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
+                                                                  eclipseState, simulator->vanguard().gridView(),
+                                                                  simulator->vanguard().cartesianMapper(), 10.0);
 
     const auto& pressures = comp.press();
     BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
@@ -590,7 +594,9 @@ BOOST_AUTO_TEST_CASE(DeckWithCapillaryOverlap)
     Opm::GridManager gm(eclipseState.getInputGrid());
     const UnstructuredGrid& grid = *(gm.c_grid());
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
+    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
+                                                                  eclipseState, simulator->vanguard().gridView(),
+                                                                  simulator->vanguard().cartesianMapper(), 9.80665);
     const auto& pressures = comp.press();
     BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
     BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@@ -649,7 +655,9 @@ BOOST_AUTO_TEST_CASE(DeckWithLiveOil)
     const UnstructuredGrid& grid = *(gm.c_grid());
 
     // Initialize the fluid system
-    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
+    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
+                                                                  eclipseState, simulator->vanguard().gridView(),
+                                                                  simulator->vanguard().cartesianMapper(), 9.80665);
     const auto& pressures = comp.press();
     BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
     BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@@ -724,7 +732,9 @@ BOOST_AUTO_TEST_CASE(DeckWithLiveGas)
     Opm::GridManager gm(eclipseState.getInputGrid());
     const UnstructuredGrid& grid = *(gm.c_grid());
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
+    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
+                                                                  eclipseState, simulator->vanguard().gridView(),
+                                                                  simulator->vanguard().cartesianMapper(), 9.80665);
     const auto& pressures = comp.press();
     BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
     BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@@ -802,7 +812,9 @@ BOOST_AUTO_TEST_CASE(DeckWithRSVDAndRVVD)
     Opm::GridManager gm(eclipseState.getInputGrid());
     const UnstructuredGrid& grid = *(gm.c_grid());
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
+    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
+                                                                  eclipseState, simulator->vanguard().gridView(),
+                                                                  simulator->vanguard().cartesianMapper(), 9.80665);
     const auto& pressures = comp.press();
     BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
     BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
@@ -900,7 +912,9 @@ BOOST_AUTO_TEST_CASE(DeckWithPBVDAndPDVD)
     Opm::GridManager gm(eclipseState.getInputGrid());
     const UnstructuredGrid& grid = *(gm.c_grid());
 
-    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(), eclipseState, simulator->vanguard().gridView(), 9.80665);
+    Opm::EQUIL::DeckDependent::InitialStateComputer<TypeTag> comp(*simulator->problem().materialLawManager(),
+                                                                  eclipseState, simulator->vanguard().gridView(),
+                                                                  simulator->vanguard().cartesianMapper(), 9.80665);
     const auto& pressures = comp.press();
     BOOST_REQUIRE_EQUAL(pressures.size(), 3U);
     BOOST_REQUIRE_EQUAL(int(pressures[0].size()), grid.number_of_cells);
diff --git a/tests/update_reference_data.sh b/tests/update_reference_data.sh
index c4b929409..87a285cbd 100755
--- a/tests/update_reference_data.sh
+++ b/tests/update_reference_data.sh
@@ -63,6 +63,8 @@ tests[spe1_brine]="flow spe1_brine SPE1CASE1_BRINE"
 tests[spe1_water]="flow_onephase spe1 SPE1CASE1_WATER"
 tests[ctaquifer_2d_oilwater]="flow aquifer-oilwater 2D_OW_CTAQUIFER"
 tests[fetkovich_2d]="flow aquifer-fetkovich 2D_FETKOVICHAQUIFER"
+tests[numerical_aquifer_3d_1aqu]="flow aquifer-num 3D_1AQU_3CELLS"
+tests[numerical_aquifer_3d_2aqu]="flow aquifer-num 3D_2AQU_NUM"
 tests[msw_2d_h]="flow msw_2d_h 2D_H__"
 tests[msw_3d_hfa]="flow msw_3d_hfa 3D_MSW"
 tests[polymer_oilwater]="flow polymer_oilwater 2D_OILWATER_POLYMER"