Heterogeneous Aquifers Research Articles

Denitrification in heterogeneous aquifers: Relevance of spatial variability and performance of homogenized parameters

Nitrate contamination of aquifers is threatening the sustainability of water resource production around the globe. Denitrification is an ubiquitous process in groundwater systems and needs to be represented in managerial and predictive efforts. Yet, fully understanding the temporal and spatial dynamics of all biogeochemical and physical properties controlling nitrate reduction at the basin scale remains a challenging task. Conceptual simplifications such as homogenization are used, but their implications on model performance are often not well assessed. This paper (1) analyzes the importance of denitrification spatial variability on the fate of non-point source nitrate contamination, and (2) assesses the performance of upscaled denitrification rate coefficient values. Heterogeneity in the hydraulic conductivity of a typical unconsolidated alluvial aquifer system and its associated uncertainty is systematically represented in our 3D numerical simulations of flow and transport. The spatial variability in the nitrate reduction capacities of the aquifer is accounted for by 3 different conceptual models: a positive and a negative correlation with the aquifer hydraulic properties and a uniform random distribution. Our results indicate that heterogeneity in the denitrification rate coefficient – as well as the model describing spatial variability – has a significant impact on concentration statistics observed at a series of extraction wells. However, we also show that, under the modeling setting used in this study (e.g., stationarity of reactivity), upscaled values of denitrification rate coefficients allow a virtually perfect reproduction of nitrate concentration statistics, a promising outcome for employing parameter estimation (calibration) of effective rate coefficients of denitrification: A detailed description of the spatial variability of the nitrate reduction capacity of the alluvial aquifer is not necessary. Yet, these homogeneous, upscaled denitrification rate coefficients result from a complex interaction between the biochemistry and the physics of the aquifer, which makes their direct estimation a challenging task.

A Dual Heterogeneous Domain Model for Upscaling Anomalous Transport With Multi‐Peaks in Heterogeneous Aquifers

AbstractNatural aquifers are often characterized by multiscale heterogeneity and complex flow networks, challenging the reliable simulation and prediction of contaminant transport. Classic stochastic transport models usually assume independent, identically distributed probability density functions when upscaling solute dynamics, but this assumption is not valid for media with multiscale heterogeneity. To address this issue, a dual heterogeneous domain model (DHDM) was proposed to quantify solute transport in heterogeneous aquifers where solute breakthrough curves (BTCs) exhibit multiple peaks and transient tailing behaviors. To efficiently solve the resultant DHDM, a Lagrangian solver was developed by combining the renewal‐reward process (to capture solute dynamics in each domain) and the state transition probability (to capture solute particle transfer between the two heterogeneous domains). The DHDM and its Lagrangian solver were then applied to simulate solute transport in three different aquifer settings (alluvial, fractured, and karst). Analyses showed that in the alluvial aquifer, early arrivals and bi‐peaks of the observed solute BTC might be caused by interconnected preferential flow paths consisting of high conductivity sediments, while the front edge or long tails of local peaks might be a result of small‐scale heterogeneity. Fracture networks and karst conduits also resulted in anomalous transport with multiple peaks, and mass transfer between the preferential flow channels and matrix led to solute retention. All these complex, coexisting anomalous transport characteristics can be simultaneously and accurately captured by the unifying DHDM approach, significantly expanding the capability of nonlocal transport models.

Seawater Intrusion in Extremely Heterogeneous Laboratory-Scale Aquifer: Steady-State Results

This work used experimental and numerical methods to investigate seawater intrusion (SWI) in a complex heterogeneous laboratory-scale aquifer. We started the analysis with a homogeneous isotropic aquifer as a reference case, then moved to heterogeneous layered aquifers. The study also investigated block-wise synthetic aquifers with different configurations. The seawater wedge toe length generally decreased under heterogeneous conditions, while the freshwater–saltwater dispersion/mixing zone generally increased when compared to the homogenous case. The saltwater–freshwater interface shows a distinct gradient change across boundaries at differing hydraulic conductivities. This was attributed to streamline refraction, which caused a reduction to the angle of intrusion when transitioning from high to low permeability zones and vice versa. The refraction also affected the mixing zone, where additional spreading was also observed when transitioning from high to low permeability zones and vice versa. When low permeability zones predominated the shoreline at the saline water boundary, this produced a shorter saline wedge in the horizontal direction, but it was more expanded vertically. This study provides insight into the general processes of SWI in heterogeneous aquifers and could be used as a basis for defining conceptual models of real-world systems. It highlights the capabilities of the image analysis to capture small perturbations.

Localization of redox-hotspots using a comparative tracer tomography approach

The identification of reactive hotspots in shallow porous aquifers is challenging, and normally requires extensive geochemical data collection combined with computationally intensive reactive transport modeling. In this study we propose an alternative, practical approach, combining conservative and reactive tracers using a tomographic setup to localize such zones.Tracer tomography combines multiple tracer tests to capture and reconstruct aquifer heterogeneities. These methods generally use “ideal” tracers that are dominantly transported by conservative advection to provide a hydraulically accurate aquifer reconstruction. We compare these conservative tracer tomography results to those obtained by using reactive tracers in a similar tomographic setup. Assuming that the conservative tracer is ideal, the differences between the two tracers will be mainly related to reactive transport. In natural systems, the biogeochemical reactions happen in localized zones, these will be the parts of the reconstructions with the greatest differences.We present the viability of the approach in a simple 2-D example, and on an aquifer analog model-based complex example. With extensive sensitivity analysis, we test the effectiveness of the method under different reaction intensities to define the chemical requirements for an ideal reactive tracer for such investigations.

Effect of permeability heterogeneity on the dissolution process during carbon dioxide sequestration in saline aquifers: two-and three-dimensional structures

The convection–diffusion process of carbon dioxide (CO2) dissolution in a saline reservoir is investigated to shed light on the effects of the permeability heterogeneity. Using sequential Gaussian simulation method, random permeability fields in two and three-dimension (2D and 3D) structures are generated. Quantitative (average amount of the dissolved CO2 and dissolution flux) and qualitative (pattern of the dissolved CO2 and velocity streamlines) measurements are used to investigate the results. A 3D structure shows a slightly higher dissolution flux than a 2D structure in the homogeneous condition. Results in the random permeability fields in 2D indicates an increase in the standard deviation of the permeability nodes enhances the dissolution efficiency, fluctuations in CO2 dissolution flux, separation between the different realizations from the same input parameters, and tendency toward more jagged convective fingers’ shape. Furthermore, the distance between the permeability nodes increases the convective fingers’ dissolution efficiency and jagged structure. The degree of freedom in 3D structures results in a higher chance of escaping from the low permeability zones and reduces the interactions between convective fingers in 3D systems. With the same variance and correlation length between permeability nodes, connectivity between high permeable zones in 3D cases are less than that of 2D cases; therefore, 2D realizations overestimate the dissolution flux of real heterogeneous 3D structures, which should be considered carefully.Article HighlightsCO2 sequestration in two and three dimensional heterogeneous saline aquifers are investigated.3D structures in homogeneous conditions show higher dissolution than 2D structures.2D realizations overestimates the dissolution flux over real heterogeneous 3D reservoirs.

Mapping discrete fracture networks using inversion of hydraulic tomography data with convolutional neural network: SegNet-Fracture

In this paper, we propose a new method to map the fracture network structure in a heterogeneous aquifer from inversion hydraulic head data measured during pumping tests in hydraulic tomography mode. This inversion tool is based on the new concept of convolutional neural networks, which provides a direct approximation to the inverse function linking fracture geometry to hydraulic data. In order to handle the highly nonlinear inverse function more effectively, an advanced neural network is developed from SegNet architecture with encoder-decoder structure, which excels in image processing to translate the water level image associated with the pumping tests at the input into a fracture map at the output. The network is trained with a synthetic dataset where the fracture structure and matrix heterogeneity are randomly generated, and the hydraulic head are obtained by solving the groundwater flow equation. The trained network accurately maps different complexity levels of fractures embedded in a matrix with heterogeneous transmissivity.As a data-driven approach, the accuracy of the mapping depends on the quantity, quality, and relevance of the synthetic dataset used in the training phase. While generating data to train the network requires effort, the trained network performs each inversion instantly. The inversion result appears to be stable even in the presence of data noise, reliably interprets the hydraulic data if they share comparable fracture and matrix properties as specified in the training models.

Integration of deep learning and information theory for designing monitoring networks in heterogeneous aquifer systems

Integration of deep learning and information theory for designing monitoring networks in heterogeneous aquifer systems

Integration of deep learning and information theory for designing monitoring networks in heterogeneous aquifer systems

Integration of deep learning and information theory for designing monitoring networks in heterogeneous aquifer systems

Aquifer flow unit analysis using stratigraphic modified Lorenz plot: a case study of Edem, eastern Nigeria

The present study aimed at evaluating aquifer flow unit of Edem by means of stratigraphic modified Lorenz plot (SMLP) employing vertical electrical sounding. The results from the vertical electrical sounding field data analysis shows that resistivity and thickness of the aquifer layer ranges from 34.8 to 67561.2 Ωm and 29.7 to 147.6 m respectively having an average aquifer thickness value of 70.248 m. Aquifer layer permeability and fractional porosity range from 1.5942 to 1657444 mD and 0.2558 to 0.3265 respectively. Results from stratigraphic modified Lorenz plot delineated a total of eight flow-units, FU1, FU2, FU3, FU4, FU5, FU6, FU7, and FU8 for the aquifer units and three different process speeds. Aquifer flow units FU1, FU3, and FU8 fall within the speed zone, flow units FU5 and FU7 are associated with the baffle zone, while flow units FU2, FU4, and FU6 make up the barrier zone of the process speeds. The aquifer flow unit speed ranges from 265406.4 to 5076424 mD with an average value of 826310.2 mD. The contour maps generated revealed the western part of the study area as characterized with high values of aquifer storage capacity and aquifer flow capacity with the highest storage and flow capacities delineated along the northwestern part of the study area. The Dykstra-Parsons coefficient of 0.99 revealed an extremely heterogeneous aquifer in the study area. This flow unit analysis will serve as a guide in accurate design and management of the aquifers.

Influence of layered heterogeneity on nitrate enrichment induced by cut-off walls in coastal aquifers

Cut-off walls have proven to be an effective countermeasure for mitigating seawater intrusion (SI). However, construction of cut-off walls may lead to upstream accumulation of nitrate contaminants in coastal aquifers. While stratified geological heterogeneities are widely present in coastal areas, how they may influence nitrate enrichment induced by cut-off walls remains poorly understood. In the present study, we employ variable-density coupling with reactive transport numerical models to understand and quantify the dynamic changes of salinity and NO3− in layered heterogeneous coastal aquifers in the presence of cut-off walls. Our comprehensive simulation results indicate that aquifer heterogeneities are a primary factor in regulating nitrate accumulation. For example, the presence of a middle low-permeability layer between two high-permeability layers reduces the influence of cut-off wall height on SI if the bottom of the wall is placed inside the low-permeability layer. While the thickness of the middle low-permeability layer has a minor impact on SI, a thicker middle low-permeability layer leads significantly aggravated nitrate contamination. Moreover, our simulations reveal that the degree of nitrate accumulation changes non-monotonically with changes in the hydraulic conductivity of low-permeability layer—the total nitrate mass (TNM) first increases and then decreases as the hydraulic conductivity of the low-permeability layer increases. We also found that the impact of heterogeneities is a function of the location and depth of the cut-off wall. With the increases of the distance between wall and the ocean, the nitrate enrichment degree of nitrate pollution in the upstream aquifer decreases while seawater intrudes further in layered aquifer (i.e. case HLH). TNM increases as the wall depth increases, but the growth rates of TNM becomes lower when the bottom of the wall is placed inside the low-permeability layer. These quantitative analyses have critical implications for the planning and design of cut-off walls to effectively control SI while at the same mitigating nitrate pollution.

Linkages between land-use change and groundwater management foster long-term resilience of water supply in California

Study RegionWe created a 270-m coupled model of land-use and groundwater conditions, LUCAS-W[ater], for California’s Central Coast. This groundwater-dependent region is undergoing a dramatic reorganization of groundwater management under California’s 2014 Sustainable Groundwater Management Act (SGMA). Study FocusUnderstanding land-use and land-cover change supports long-term sustainable water management. Anthropogenic water demand has depleted groundwater aquifers worldwide, while future water shortages will likely affect land-use change, creating system feedbacks. Our novel participatory approach fused changes in land-use and associated water use from county-scale data to local water agencies’ estimates of total sustainable supply, scaling up local hydro-geologic knowledge from heterogeneous aquifers and diverse management approaches to a regional level. We assessed five stakeholder-driven scenarios with the same historic rates of urban and agricultural land-use change, but different water and land-use management, analyzing how management strategies altered both the spatial pattern of development and subsequent water sustainability from 2001 to 2061. New Hydrological Insights for the RegionTransformative strategies using demand-side interventions that coupled water availability to land-use more effectively achieved long-term sustainability than adaptive strategies using supply-side interventions to increase water supplies. Limiting water withdrawals within SGMA regulated basins resulted in leakage of development into unregulated basins, increasing groundwater pumping there. Protecting ecosystems, farmlands, and recharge areas from development reduced leakage into undeveloped basins without negatively affecting water sustainability.

Seawater Intrusion-Retreat Processes and Groundwater Evolution in Intruded Coastal Aquifers with Land Reclamation: A Case Study of Eastern Jiangsu, China

Abstract The hydrogeochemical analysis and numerical simulation were conducted to explore the impact of land reclamation on submarine groundwater discharge and groundwater evolution in the seawater intrusion-retreat processes in eastern Jiangsu, China. The hydrogeochemical results indicated that the origin of salt in most groundwater is saltwater mixing, and weathering-dissolution of minerals enriches the salt in fresh groundwater. The evolution path of coastal groundwater in Jiangsu is Ca-Na-HCO3-Cl → Ca-Na-Cl-HCO3 → Na-Ca-Cl-HCO3 → Na-Mg-HCO3 → Na-Mg-HCO3-Cl → Na-Cl-HCO3 → Na-Cl, and the difference in groundwater evolution in different regions is obvious. The numerical simulation revealed that the land reclamation has blocked groundwater exchange and restrained the seawater intrusion. The sensitive analysis indicated that the permeability and the slope of the land reclamation have dominated the salinity distribution. And the aquifer heterogeneity promotes groundwater evolution and leads to the irregular salinity distribution of the groundwater. Moreover, the groundwater evolution in the shallow aquifer with the land reclamation shows that the residual seawater is gradually desalinating and discharged into the sea. The land reclamation is salinized by saltwater discharge and seawater intrusion together and then subsides into a stable saltwater wedge. The groundwater salinity distribution in Jiangsu is not fit seawater intrusion, which is caused by the coastline changes due to anthropogenic land reclamation and natural activities.

New challenges in the study of coastal aquifers from numerical modeling. The case of the Motril-Salobreña aquifer

Coastal aquifers are frequently complex systems with highly heterogeneous geological characteristics, specific flow patterns, presence of fluids with different densities, high risk of contamination both by salinization and by other pollutants, and highly modified by anthropic activity. Given this situation, numerical modeling becomes the most appropriate tool to determine the potential impact of global change on marine intrusion in this kind of aquifers. The mathematical models traditionally used in coastal aquifers have been those of flow and mass transport with variable density, which allow obtaining a distribution of salinities in the aquifer and reproducing the flow pattern in the area of discharge to the sea. In addition to these models, another type of numerical modeling can be applied that could also provide information on specific aspects of this type of aquifers. The study of the Motril-Salobreña coastal aquifer shows, for example, how heat transport models allow quantifying the recharge that occurs from rivers that frequently interact with coastal aquifers; as well as the age transport models provide data on the preferential groundwater circulation paths and its residence time, very relevant information in these frequently heterogeneous and anisotropic aquifers. Future challenges are directed towards determining the effects of sea level rise on marine intrusion and establishing the degree of equilibrium of the freshwater-saltwater contact with the current situation.

Hydro-thermo-chemo-mechanical modeling of carbon dioxide injection in fluvial heterogeneous aquifers

Understanding geological sequestration of carbon dioxide (CO2) requires fully-coupled simulation of thermal-hydrological-mechanical-chemical (THMC) processes within well-characterized subsurface heterogeneity. A CO2 injection pilot project was performed in the Cranfield site, Mississippi, USA, where subsurface heterogeneity reflects fluvial deposition. During the project, CO2 was injected through an injection well and was monitored using two observation wells. We incorporate high-resolution three-dimensional heterogeneity model of the site into multiphase and multi-component flow and transport models. We validate our models and evaluate how bottom-hole pressure (BHP) in injection well, CO2 breakthrough times at observation wells, and efficiency of trapping mechanisms (i.e., dissolution, snap-off trapping, mobile CO2) are controlled by (1) non-isothermal CO2 injection, (2) geochemical reactions, (3) capillary pressure heterogeneity, (4) geomechanical effects, and (5) permeability enhancement close to the injection well. Results suggest that neglecting thermal effects lowers BHP, shortens breakthrough times, overestimates dissolution, and underestimates snap-off trapping. The BHP remains unchanged, breakthrough times are overestimated, and solubility and snap-off trapping are underestimated when geochemical reactions are ignored. Ignoring capillary pressure heterogeneity results in underestimation of BHP and dissolution and overestimation of breakthrough times. Ignoring geomechanical process results in lower BHP, shorter breakthrough times, increased dissolution, and decreased snap-off trapping. Ignoring permeability enhancement results in underestimation of breakthrough times and solubility trapping, and overestimation of snap-off trapping. Our results have important practical implication for designing effective field scale experiments and numerical simulations of geological carbon sequestration. Other multiphase flow problems in which subsurface heterogeneities are the main controlling factor (e.g., hydrogen storage) may benefit from this work.

Exploration of Ordovician limestone aquifer heterogeneity with tomographic water releasing tests

Water inrush is a typical geological disaster in the process of underground mining. It is crucial to depict potential water inrush areas for mining safety. In this study, mining activity in Dongtan coal mine is under the water inrush threat from deep Ordovician limestone confined aquifers. Tomographic water releasing tests and simultaneous successive linear estimator (SimSLE) were conducted to characterize aquifer heterogeneity by stimulating hundreds of square kilometers. An artificial injection well was firstly introduced to alleviate the influence from identified recharge channels. The estimated parameter distributions indicated that the study area was divided into two separately hydrogeological compartments. Transmissivity (T) and storativity (S) in the north were relatively higher due to water-conducting structures, and more attention should be paid to this region which may develop into water inrush incidents. The southern area is occupied with lower T and S, which indicated relatively worse flow conditions. Auxiliary drawdown recovery data, borehole chemistry variations, drilling fluid losses, and existing geological features further confirmed the estimated parameter distributions. In the end, we concluded that Hydraulic Tomography technique (HT) could obtain aquifer parameter information in large scales and provide supports to formulate reasonably effective water prevention and control schemes before mining deep coal seams.

Entity aware sequence to sequence learning using LSTMs for estimation of groundwater contamination release history and transport parameters

Groundwater contaminant sources identification and parameter estimation using the simulation–optimization (S/O) approach require numerous runs of the computationally expensive simulation model through the optimization algorithm. The computational cost can be effectively reduced by using a surrogate model which can accurately approximate the simulation model. With the advent of deep learning, Long Short-Term Memory (LSTM) networks, which are suitable to learn sequential data, are being increasingly applied to regression problems involving time dependencies. However, for the simultaneous contaminant source identification and parameter estimation problem, the surrogate model requires to establish a relationship between release histories at the source locations to concentration measurements at observation points subject to given values of aquifer parameters. In this study, a novel deep neural network framework Entity Aware Sequence to sequence learning using Long Short-Term Memory (EAS-LSTM) is proposed as a surrogate model which takes both sequential i.e., release histories at the source locations at different stress periods and static variables i.e., transport parameters as inputs to predict breakthrough curves (BTCs) at observation points. The proposed surrogate model is applied to a heterogeneous field scale aquifer with 4 zones. The Mean Squared Error (MSE) using EAS-LSTM is 0.0033 ppm2 which is significantly better than that of Kriging and Support Vector Regression (SVR) which are 0.048 ppm2 and 1.302 ppm2 respectively. Comparison of EAS-LSTM model performance with Kriging and SVR based models demonstrates its higher accuracy. Further, the optimization algorithms for inverse modelling, the performances of Multiverse Optimizer (MVO), Grey wolf optimization (GWO) and Particle swarm optimization (PSO) are investigated. It is observed that the combination of EAS-LSTM and MVO provides better results in comparison to other surrogate simulation optimization (SSO) models.

Aquifer Characterization Using Fiber Bragg Grating Multi-Level Monitoring System.

Groundwater responses measured from multiple wells at different depths are essential for delineating the aquifer heterogeneity using hydraulic tomography (HT). In general, conducting HT requires many wells because traditional well monitoring is usually partially open at a specific depth interval or is fully penetrating. Accordingly, conducting an HT survey is typically costly and time-consuming. To tackle these issues, a new multi-level monitoring system (MLMS) for the HT survey was developed using the fiber Bragg grating (FBG) technique. This FBG MLMS could collect the depth-discrete groundwater observations from a fully penetrated 2-inch well. Three field campaigns were conducted to validate the capability of the FBG MLMS for HT surveys. The results show that the accuracy and stability of this MLMS are reliable and that FBG MLMS is beneficial for conducting an HT survey. Specifically, compared to the traditional monitoring well in an injection event, this FBG MLMS can concurrently cause an increase in the number of cross-hole tests several times and collect many more head observations than the standard methods, resulting in the observed flow fields efficiently reaching ergodic conditions and effectively improving the accuracy of the estimated hydraulic heterogeneity. Therefore, the FBG MLMS could be an alternative MLMS for efficiently and economically conducting an HT survey.

Mechanistic modeling and measurement of foamed gas flow resistance in fractures

Geologic heterogeneity is one of the important factors controlling migration and distribution of nonaqueous phase liquids in the subsurface. Heterogeneity, such as fractures, complicates the effectiveness of environmental remediation, and uncharacterized aquifer heterogeneity causes significant uncertainties in the removal of nonaqueous phase liquids. Foaming an injection gas is an effective method to alleviate the challenges in accessing heterogeneous media during remediation processes because foam is nonNewtonian and redirects gas and liquid flow. Foamed gas may also assist in increasing the storage of carbon dioxide in saline formations. This paper investigates foam flow behavior in microscale fractures and develops a mechanistic transient foam flow model using the population balance method. Microscale experiments in fractures with apertures of 25 and 88 μm are used to validate the model for pressure drop, gas saturation, and bubble texture. Key differences related to modeling foam in fractures are the potential for continuously varying gas–liquid curvature in fractures and the relationship of this curvature to apparent foam viscosity. Incorporation of a local foam flow resistance factor is important to representing flow physics accurately. Hence, the results of this study provide foundational understanding for upscaling foam rheology in different sized fractures to ensure successful operations in soil/aquifer remediation and CO2 storage in saline formations.

Deep learning based optimization under uncertainty for surfactant-enhanced DNAPL remediation in highly heterogeneous aquifers

The optimization for surfactant-enhanced aquifer remediation (SEAR) of dense non-aqueous phase liquid (DNAPL)-contaminated aquifers is usually accompanied by uncertainties, which may arise from the characterization of complex aquifer heterogeneity and DNAPL source zone architecture (SZA) due to measurement sparsity. Optimization under uncertainty is computationally expensive as it involves an enormous number of model runs. The huge computational burden can be alleviated by utilizing a surrogate model for repeated model evaluations. However, most of the developed surrogates are often limited to low-dimensional optimization problems that only consider simplified aquifer heterogeneity. In this study, we developed a multi-objective simulation-optimization framework to optimize the SEAR schemes considering the characterization uncertainties from both highly heterogeneous aquifer permeability and complex SZA. A fast-to-run convolutional neural network (CNN)-based surrogate model was developed to approximate the high-dimensional and highly complex input-output mapping of the DNAPL multiphase flow simulation model. We first used the rejection sampling strategy to generate random realizations of permeability and SZA conditioning on their limited measurements and then formulated a multi-objective optimization under uncertainty problem based on these realizations. The developed 3-D CNN was trained and used as the surrogate for repeated model runs in optimization to identify the optimal SEAR schemes under uncertainty. A 3-D numerical experiment was used to test the performance of the CNN-based simulation-optimization framework. Comprehensive analysis on the obtained Pareto fronts demonstrates that the proposed framework can efficiently identify reliable Pareto-optimal solutions with a 99.8% speedup compared to the traditional optimization coupled with the forward model. Moreover, the optimization considering multiple realizations enables us to perform the risk assessment to locate the risk zone where the NAPL phase possibly exists after remediation, which provides useful information for decision-making.

A Sequentially- Hybridized Locally Conservative Non-conforming Finite Element Scheme for Two-phase Flow Simulation through Heterogeneous Porous Media

Accurate two-phase flow modeling has a crucial role in detecting pollution in aquifers and/or in the simulating immiscible flow through porous media. In this research, a new hybrid numerical method based on the lower order non-conforming finite element method (NCFEM) and interior penalty discontinuous Galerkin (DG) method is proposed to simulate two-phase incompressible flow in heterogeneous porous media. The pressure and transport saturation equations are discretized by using the non-conforming Crouzeix-Raviart (CR) finite element and the symmetric weighed interior penalty discontinuous Galerkin (SWIPM). As a result of determining the degrees of freedom at the midpoint of neighboring element edges, a consistent set of pressure and velocity fields is obtained via the NCFEM. An H(div) projection based on the Raviart-Thomas (RT) element is implemented to improve the results resolution and preserve the continuity of the normal component of the velocity field. Besides, the spurious numerical oscillations of the saturation solution at the end of each time step are removed by using a novel vertex-based slope limiter, named the modified Chavent-Jaffre limiter. The accuracy and performance of the proposed numerical method are assessed by solving the Buckley-Leverett and McWhorter-Sunada benchmark problems, along with two test problems associated with highly heterogeneous porous media. The numerical results indicate that the proposed scheme has a remarkable potential to accurately capture the shock front and sharp interface zone of immiscible phases in the heterogeneous aquifers and other porous media.