Numerical modeling of reservoir stimulation using matrix acidizing requires a thorough understanding of chemical reactions, fluid flow, heat transfer, and their interactions in porous rock. This research introduces a robust in-house mesh-free code utilizing Element Free Galerkin (EFG) method to address Chemical, Hydro-Chemical (HC), and Thermo-Hydro-Chemical (THC) problems, particularly used for simulation of matrix acidizing operations with a two-scale continuum approach. Incorporating thermal interactions improves the realism of simulations by considering the heat from exothermic reactions. This study also introduces a combination of Equivalent Continuum Method (ECM) with EFG method to examine the impact of the rock microcracks on the matrix acidizing process. The verification of the developed EFG code is carried out first by solving an analytical chemical transport problem (C). Then two HC, and one THC matrix acidizing problems are solved and the obtained results confirm its ability to model acid flow, chemical and thermal effects, highlighting its efficacy for simulating various matrix acidizing problems. In the EFG-ECM simulation of the effects of microcracks in a coupled THC study, 0.15%, 0.71% and 12.28% reductions in Pore Volume to Break Through (PVBT) were observed for 0o, 45o and 90° crack orientations respectively, in comparison with similar cases without thermal effects. The results highlight a minimal temperature impact on the PVBT for 0o and 45° crack inclination cases, due to unchanged dissolution pattern; however, temperature variation influences acid kinetics and dissolution pattern in the vertical cracks case, resulting in a noticeably lower PVBT.
Salinity-induced density variations significantly impact flow dynamics and the migration of petroleum hydrocarbons (PHC). Although density-dependent flow and PHC transport are extensively modeled independently, their coupled interaction remains insufficiently understood, particularly under variably saturated conditions. In this study, a coupled density-dependent flow and multi-component transport model is developed to simulate the impact of salinity on BTEX (benzene, toluene, ethylbenzene and xylene) migration. The influence of density-dependent flow is analyzed across various soil textures, salt concentrations and source strengths under variable saturated conditions. Vertical flow is simulated in one-dimensional domain, while transport equations are solved in two dimensions. The equations are discretized using fully-implicit finite difference scheme with nonlinearities addressed using Picard iteration. The discretized equations for flow and transport are solved using Tridiagonal Matrix Algorithm (TDMA) and Bi-Conjugate Gradient Stabilized method (BiCGSTAB). Model results show that density variation significantly affects flow dynamics and BTEX transport in sand and loam and has negligible effect in clay. At an inlet salt concentration of 300 g/L, density variations increase Darcy flux by approximately 19% in sand and 17% in loam after 40 days. Density effects increase BTEX dissolution rate and enhance the migration depth of BTEX by 10 - 40% in sand. The impact of salinity on BTEX transport intensified with the increase in the BTEX source zone strength. BTEX transport is found to be sensitive to soil hydraulic parameters and biodegradation kinetics. These findings highlight the significance of considering the density-dependent flow while predicting BTEX transport under salt co-transport scenarios.
Quantifying pore-scale fluid displacement mechanisms in CO2/brine systems is essential for predicting multiphase flow and trapping during CO2 storage. In this study, we imaged the steady-state flow of brine and CO2 in a water-wet Bentheimer sandstone sample. We devised an experimental method based on differential imaging to examine the pore-scale behaviour of CO2 within the system at reservoir conditions during drainage. This allowed us to measure relative permeability and capillary pressure while observing the evolution of the gas ganglia in the pore space at different fractional flows under capillary-dominated conditions.The measured CO2 relative permeability remained low during the initial stages of drainage across a broad saturation range and only increased modestly to 0.24 at 100% CO2 injection, with a gas saturation of 0.57. Initially, CO2 appeared as small and disconnected ganglia occupying the largest pores and progressively invaded smaller pores and throats with increasing gas fractional flow, leading to partial connectivity.The interface exhibited predominantly positive mean curvature, consistent with a water-wet system. Capillary pressure estimated from interfacial curvature showed good agreement with independent porous-plate measurements reported in the literature. Overall, the connectivity of the CO2 phase, quantified by Euler characteristic, remained limited during the displacement process contributing to the measured low relative permeability. These observations are consistent with a capillary-dominated displacement process with local trapping by Roof snap-off.
This paper presents a new model for analysing the advection-dispersion of a contaminant in a pipeline under unsteady flow conditions. Compared to other papers in the scientific literature, the analysis focuses on fast transient, rather than steady, quasi-steady or slow transient conditions. The new model uses the method of characteristics (MOC) and the Crank-Nicolson finite difference scheme, to couple flow field with contaminant advection and to model dispersion, respectively. For fast transient analysis, the developed MOC-based scheme accounts for the unsteady friction effect through a simplified acceleration based approach. Firstly, the model is tested against advection and dispersion under steady flow while compared to other numerical schemes. The model’s performance is evaluated in terms of accuracy, computation time and stability. Comparisons with benchmark models and exact solutions prove the effectiveness of the present numerical coupling methodology at preserving the shape of the concentration profile (more than 97%) and its reduced proneness to instability. The model is applied to a challenging test case under fast transient flow conditions. Therefore, numerical diffusion as well as changes in contamination profiles induced by the variation of the flow velocity is emphasized and thoroughly discussed.
Injecting hot flue gas into coal seams is a highly promising technology that can enhance the permeability of coal seams and achieve carbon sequestration. When hot flue gas is injected into water-containing coal seams or when it carries water vapor, the solution within the reservoir forms a weakly acidic environment, causing the minerals in the coal to dissolve, thereby altering the pore structure. However, observing the structural dynamic changes of coal samples during the dissolution process through experiments is usually costly and is also limited by time resolution and duration. To address this issue, we developed a phase field model that combines heat, chemical reactions, and diffusion effects to study the dissolution process in mineral-containing porous media. The model results were verified against the data obtained from static coal experiments conducted under specific temperature and pressure conditions in saturated carbon dioxide solutions. The structural data of the coal samples were obtained through X-ray computed tomography (CT) technology. Using the original slice structure of the coal as the model geometry, the entire dissolution process was simulated. The results show: (1) The reaction rate-time curve roughly follows a negative exponential decay pattern, and there are obvious turning points, presenting a nearly stepped curve characteristic; (2) Three different dissolution modes were determined: fully fracture-dominated type; fracture-dominated and pore secondary type; and pore-dominated and fracture secondary type. These findings provide theoretical insights into the mineral dissolution patterns at the fracture scale in complex porous media.
Stratification in river–reservoir systems receiving multiple inflows emerges from the interaction between buoyancy-driven density gradients and shear-induced turbulent mixing. In such systems, asymmetric inflows and hydraulic regulation can generate pronounced lateral heterogeneity, which may limit the representativeness of conventional cross-sectionally averaged stability assessments.This study investigates stratification regimes in a multi-inflow river–reservoir system using high-resolution in-situ hydrodynamic and thermal observations. A sectional bulk Richardson number approach was applied to resolve transverse variability in shear–buoyancy balance across multiple transects. By evaluating stability independently at left, center, and right sections, the method identifies localized stratification structures that may not be fully captured by cross-sectional averaging.Three stratification regimes were observed: shear-dominated mixed conditions, transitional states, and buoyancy-dominated stratified conditions. Periods of strong density forcing exhibited enhanced lateral variability in stability, indicating that inflow asymmetry and hydraulic regulation can influence regime development in addition to seasonal thermal forcing.
Although climate fundamentally shapes watershed evolution, our understanding of the underlying processes and their relationships with climatic factors remains limited. Recent advances in landscape evolution models now enable the quantitative analysis of fluvial system development. In this study, we employ a landscape evolution model to simulate watershed evolution, specially examining how climatic variables impact watershed geomorphic system and river network morphology. The results demonstrate that climate humidity is a primary driver of watershed geomorphology and river network morphology. Characteristic parameters-including the box-counting fractal dimension, geomorphic fractal dimension, main channel length, and drainage density-exhibit notable positive correlations with humidity, whereas the hypsometric integral is negatively correlated. The functional forms of these relationships vary: the hypsometric integral decreases linearly with humidity, while the box-counting fractal dimension and main channel length follow convex trends. Furthermore, drainage density displays concave curves under conditions of inactive tectonic uplift, transitioning to convex curves when tectonic uplift is active. The maximum stream order of the generated river networks increases with humidity; specifically, the geomorphic fractal dimension ranges from 1.2 to 1.4 for third-order streams, 1.4 to 1.7 for fourth-order streams, and 1.9 to 2.0 for fifth-order streams. Within a given stream order, both the bifurcation ratio and the length ratio positively correlate with humidity, with higher-order rivers showing a narrower range of variation. Finally, the bifurcation ratio follows an approximately linear relationship with humidity, whereas the length ratio conforms to an S-curve.
Multiscale heterogeneity in hydraulic conductivity and the sparse, finite-support nature of well measurements complicate the reconstruction of contaminant plume dynamics in natural aquifers. While data-driven reduced-order methods provide a computationally efficient alternative to fully resolved simulations, their sensitivity to observational design remains insufficiently examined. This study evaluates the ability of Dynamic Mode Decomposition (DMD) to recover and extrapolate concentration fields in heterogeneous porous formations using spatially sparse, volume-averaged observations. High-fidelity simulations serve as reference solutions across varying heterogeneity levels and Péclet numbers, enabling systematic assessment of how reconstruction accuracy depends on transport regime, sampling density, and measurement support scale. DMD captures dominant plume structures when observational density is sufficient to resolve the underlying heterogeneity. Performance is strongly governed by the interaction between advective–dispersive dynamics and observational scale: increased sampling density and larger measurement support enhance stability by attenuating small-scale variability, whereas sparse observations in strongly advective, highly heterogeneous settings reduce extrapolation reliability. These results quantify how observational density and sampling support influence data-driven plume reconstruction, offering preliminary guidance for monitoring strategy design under similar idealized transport conditions.
Large-scale transient groundwater monitoring using machine learning (ML) remains challenging, primarily due to limited observational data. Here we take the contiguous United States (CONUS) as a case study. Despite records from over one million water table depth (WTD) monitoring wells, groundwater observations remain spatially and temporally sparse. We systematically quantify multiple dimensions of these groundwater observation gaps. To evaluate the impact of data limitations on model performance, we train a random forest model for WTD seasonality estimation over the CONUS based on a balanced dataset that contains training and test sets of equal sizes across different land cover types. Overall, the model performs relatively well, with a median Spearman’s ρ of 0.53, but still fails to reproduce seasonal WTD variations in many locations, especially those impacted by human activities. Comparison with a random forest trained on physically-based simulation results at the same grid cells suggests the complexity in modeling human-impacted groundwater systems. Previous hydrometeorological conditions are identified as the most important input for all land cover types. A perturbation experiment further shows that, under the current groundwater monitoring network across the CONUS, the model predictive performance is more sensitive to additional observations in grasslands and farmlands than those in other land cover types, offering empirical guidance for future well installation. Our study provides insight into how data challenges affect ML model predictive performance and highlights the importance of optimizing and developing global groundwater monitoring networks.
Joule-Thomson (JT) cooling during CO2 injection poses significant thermal and operational challenges in depleted reservoirs. This study performs fully coupled wellbore-reservoir simulations to investigate transient thermal-hydraulic behavior during CO2 injection and to explicitly compare depleted gas and brine-filled reservoir conditions. Results show that depleted reservoirs experience pronounced JT cooling, with near-wellbore temperatures dropping to approximately 13 °C within the early injection stage. The cooling front propagates rapidly during the initial period and then transitions to a diffusion-controlled regime, with the 38 °C isotherm following a logarithmic radial expansion trend. A high-density liquid CO2 zone forms near the wellbore due to combined pressure buildup and cooling effects. In contrast, the brine-filled reservoir exhibits smaller temperature reductions and smoother pressure evolution, owing to its higher initial formation pressure. Parametric analysis indicates that injection rate strongly influences the intensity of transient cooling, with higher rates amplifying the early-stage temperature decline. Increasing the injection temperature mitigates but does not eliminate JT cooling under low reservoir pressure conditions. These findings highlight the dominant role of the thermodynamic state of the initial reservoir in governing injection-induced cooling and define operational strategies to maintain thermal integrity during CO2 storage in depleted gas reservoirs.
High-temperature geothermal and shale gas extraction rely on gas flow through artificially generated shear and tensile rock fractures. Existing experiments mostly focus on low-Re water flow in single fracture types, while systematic high-Re nitrogen gas comparisons and predictive models for granite fractures remain limited. Triaxial flow tests were performed on two shear and three tensile granite specimens under 1–40 MPa effective confining stress. 3D laser scanning quantified three roughness indices (Rp, Rrms, Ra), and all samples were strictly aligned to eliminate offset between matching fracture surfaces. A compressibility-modified Forchheimer equation and non-Darcy factor E were adopted to divide flow into laminar (E = 10%)), transitional (10%<E < 90%)), and turbulent (E = 90%)) regimes. Results show shear fractures possess larger hydraulic apertures, higher flow velocities and lower critical Reynolds number Rec, dominated by transitional and turbulent flow; tensile fractures mostly maintain laminar and transitional states without turbulence. Two independent empirical models were established in this work. First, an exponential function is developed to characterize the attenuation trend of nonlinear coefficient B' with growing hydraulic aperture, which achieves better fitting performance than conventional power-law expressions. Based on this exponential model of B', a novel predictive Rec model is further derived to reproduce the unique 'increase-then-stabilize/decrease' variation of critical Reynolds number with increasing confining stress. Second, a negative allometric model is proposed to quantify the logarithmic correlation between critical pressure gradient and hydraulic aperture. Cross validation against publicly available gas flow experimental data verifies the accuracy of all proposed models. This work clarifies hydro-mechanical disparities between two fracture types and supplies practical prediction tools for geothermal and shale gas engineering.
