Reactive transport models (RTM) are critical tools for understanding complex fluid-rock interactions in fractured rock systems. However, applying RTM in such systems remains challenging due to the wide disparity of scales spanning from millimeters to kilometers. While existing RTM effectively address fracture-matrix interactions at smaller scales, few can handle irregular fracture networks on larger scales while explicitly incorporating matrix processes. To bridge this gap, we introduce a novel discrete fracture-matrix (DFM) reactive transport model based on the MIN3P code, utilizing an anisotropic quadrilateral mesh. This approach enables coarser discretization along fractures (advection-dominated zones) and refined discretization perpendicular to fractures (diffusion-dominated zones), significantly improving computational efficiency without compromising accuracy. The model's capabilities are demonstrated through simulations of conservative tracer transport and dissolved oxygen migration in fractured crystalline rock at both intermediate (hundred-meter) and large (kilometer) scales. Comparative analyses with traditional triangular mesh methods reveal that the proposed approach achieves comparable or superior accuracy while drastically reducing computational demands. The model's ability to efficiently simulate large-scale fractured rock systems makes it a powerful tool for applications such as assessing geochemical stability in crystalline rock with extensive fracture networks.

Agricultural soils significantly contribute to greenhouse gas (GHG) emissions, influencing global climate change through natural microbial processes and human activities. Soil structure plays a crucial role in regulating these processes, especially in driving the formation of GHG hotspots and hot moments. We developed a new multi-domain reactive transport model to assess the mechanisms governing GHG emissions from macroporous agricultural soils, focusing on N2O hotspots and hot moments. The multi-domain model adequately reproduced observed soil moisture, pore gas O2 and GHG levels, and short- and long-term fluxes at an experimental site in eastern Ontario, Canada, outperforming other modeling approaches (i.e., uniform and dual porosity/permeability models). Building on the calibrated model, we focused on investigating the development of N2O hotspot and hot moments. While consistent with conventional assumptions that N2O production is concentrated in organic soils and subsoils, our results further reveal the fine-scale spatiotemporal variability of hotspots driven by localized variations in moisture, anaerobic conditions, and nutrient and organic carbon availability. N2O hot moments occurring during periods of heavy rainfall were closely linked with N2O hotspots in organic soils. A combination of sensitivity analysis and structural equation modeling suggests that the formation of N2O hotspots and hot moments is strongly influenced by physical and geochemical parameters, such as immobile matrix radius, gas saturation within the mobile matrix, and the ratio of N2O production to consumption. While aligning with field data, the model's complexity poses challenges, including non-uniqueness due to uncertain inputs, subdomain characterization, and gas/solute transfer between the subdomains. Nevertheless, the multi-domain model provides valuable insights that can guide future research into how complex soil structural characteristics influence GHG cycling and emissions from agricultural soils.

Agricultural soils significantly contribute to greenhouse gas (GHG) emissions, influencing global climate change through natural microbial processes and human activities. Soil structure plays a crucial role in regulating these processes, especially in driving the formation of GHG hotspots and hot moments. We developed a new multi‐domain reactive transport model to assess the mechanisms governing GHG emissions from macroporous agricultural soils, focusing on N 2 O hotspots and hot moments. The multi‐domain model adequately reproduced observed soil moisture, pore gas O 2 and GHG levels, and short‐ and long‐term fluxes at an experimental site in eastern Ontario, Canada, outperforming other modeling approaches (i.e., uniform and dual porosity/permeability models). Building on the calibrated model, we focused on investigating the development of N 2 O hotspot and hot moments. While consistent with conventional assumptions that N 2 O production is concentrated in organic soils and subsoils, our results further reveal the fine‐scale spatiotemporal variability of hotspots driven by localized variations in moisture, anaerobic conditions, and nutrient and organic carbon availability. N 2 O hot moments occurring during periods of heavy rainfall were closely linked with N 2 O hotspots in organic soils. A combination of sensitivity analysis and structural equation modeling suggests that the formation of N 2 O hotspots and hot moments is strongly influenced by physical and geochemical parameters, such as immobile matrix radius, gas saturation within the mobile matrix, and the ratio of N 2 O production to consumption. While aligning with field data, the model's complexity poses challenges, including non‐uniqueness due to uncertain inputs, subdomain characterization, and gas/solute transfer between the subdomains. Nevertheless, the multi‐domain model provides valuable insights that can guide future research into how complex soil structural characteristics influence GHG cycling and emissions from agricultural soils.

Enhanced Rock Weathering (ERW), in which crushed basaltic rocks are spread on croplands, has emerged as a promising carbon dioxide removal (CDR) approach to mitigate climate change impacts. Important known constraints on weathering rates include temperature, humidity, and feedstock grain size. However, the quantitative prediction and optimization of CDR is currently limited by uncertainty in the processes and rates governing weathering and export. Here, we propose to evaluate the product of effective groundwater recharge and dissolved inorganic carbon (DIC) concentrations as a measure of CDR. Since maximum DIC concentrations in pore water are controlled by soil gas PCO 2 and achievable Ca and Mg concentrations from weathering, we define this CDR flux as the “carrying capacity”. We consider the onset of precipitation of secondary Ca-carbonate minerals in solution due to the accumulation of ERW reaction products in the shallow soil pore water as the upper limit for the effective generation of CDR. Our results therefore present values of “maximum efficient CDR”, yielding CDR export to groundwater in the absence of carbonate mineral precipitation, as a function of effective groundwater recharge. Extending the carrying capacity concept to global croplands highlights the potential importance of groundwater recharge in determining regions with highest ERW potential. Given the simplifying assumptions in our assessment, we estimate a global CDR potential of 0.15 and 0.85 Gt CO 2 yr − 1 . Our results indicate that regions with high groundwater recharge and feedstocks rich in leachable Mg provide the highest potential for efficient CDR generation, assuming feedstock dissolution is not limiting. Our analysis does not account for CDR losses in the near field or far-field zone, for example due to nitrification or the release of stored acidity, but also omits soil exchange and other reactions that may restrict calcite precipitation and thus lead to higher maximum efficient CDR.

This paper presents and discusses the results obtained by the participants to the benchmark described in de Hoop et al, Comput. Geosci. (2024). The benchmark uses a model for CO2 geological storage and focuses on the coupling between two-phase flow and geochemistry. Several test cases of various levels of difficulty are proposed, both in one and two spatial dimensions. Six teams participated in the benchmark, each with their own simulation code, though not all teams attempted all the cases. The codes used by the participants are described, and the results obtained on the various test cases are compared, as well as the performance of the codes. It is shown that the results obtained are widely consistent, giving a good level of confidence in the outcome of the benchmark. The general complexity of two-phase flow coupled with chemical reactions altering porous media means that some differences between the codes remain. Besides, from the convergence study, it is clear that the two-dimensional problem has a relatively high sensitivity to a spatial resolution which adds to the complexity.

Oxygen (O 2 ) availability in soils is vital for plant growth and productivity. The transport and consumption of O 2 in the root zone is closely linked to soil moisture content, the spatial distribution of roots, as well as structure and heterogeneity of the surrounding soil. In this study, we measure three‐dimensional root system architecture and the spatiotemporal dynamics of soil moisture (θ) and O 2 concentrations in the root zone of maize ( Zea mays ) via non‐invasive imaging, and then construct and parameterize a reactive transport model based on the experimental data. The combination of three non‐invasive imaging methods allowed for a direct comparison of simulation results with observations at high spatial and temporal resolution. In three different modeling scenarios, we investigated how the results obtained for different levels of conceptual complexity in the model were able to match measured θ and O 2 concentration patterns. We found that the modeling scenario that considers heterogeneous soil structure and spatial variability of hydraulic parameters (permeability, porosity, and van Genuchten α and n ), better reproduced the measured θ and O 2 patterns relative to a simple model with a homogenous soil domain. The results from our combined imaging and modeling analysis reveal that experimental O 2 and water dynamics can be reproduced quantitatively in a reactive transport model, and that O 2 and water dynamics are best characterized when conditions unique to the specific system beyond the distribution of roots, such as soil structure and its effect on water saturation and macroscopic gas transport pathways, are considered.

Two-dimensional reactive transport models, one with a simplified root system and the other accounting for dynamically evolving root architecture, were constructed to examine the influence of model complexity on capturing the effect of soil-root dynamics relating to the Oxalate Carbonate Pathway (OCP) of the Iroko tree over 170 years. Oxidation of oxalate from fallen tree tissue by soil bacteria enables local soil pH increase, leading to the sequestration of atmospheric carbon in carbonate minerals (calcite) in the shallow soil surrounding the tree. Simulations of both root models corroborate previous one-dimensional models of the OCP focused on Ca and C mass balance, where high weathering rates of Ca-containing silicate minerals in bedrock, along with contributions from groundwater, provided sufficient Ca for precipitation of observed quantities of calcite. Both simulations demonstrate the development of a distinct high pH zone where oxalate is oxidized, Ca accumulates, and calcite precipitates (OCP zone); and a low pH zone where roots collect Ca, later returned to the top soil as calcium oxalate (Total Root Extent/TRE zone) via litterfall. While the extent of OCP zone development near the ground surface was very similar between simulations, differences in localized root water uptake between the two approaches resulted in variation in water and solute transport and influenced the geometry of the OCP zone at depth, with implications for calcite precipitation in the soil. Trends in CO2 and O2 partial pressures in the OCP zone were mirrored in the TRE zone, suggesting linkage between the two zones with regard to gas transport. Near the end of the tree's lifespan, results indicate that soil permeability decreases due to calcite precipitation may limit O2 ingress and availability in the shallow soil, while trapping CO2 released from the oxidation of organics in the shallow soil, with implications for the long-term sustainability of the OCP itself.

Agricultural soil is the domain where biogeochemical C-N turnover is affected by atmospheric and hydrologic processes, playing a vital role in anthropogenic greenhouse gas (GHG) emissions and thus affecting global climate change. GHG cycling and emissions are further complicated by the characteristics of soil structure, in particular macropores (i.e., preferential flow pathway). In this study, the processes of non-equilibrium gas diffusion and gas exchange were implemented into an existing dual-permeability flow and reactive solute transport model. The model was then used in a sensitivity analysis to evaluate the suitability of the modeling approach to account for development of anaerobic hotspots, with an emphasis on assessing N2O production and emissions. The sensitivity analysis combined with characteristic time scale analysis suggests that the development of anaerobic hotspots is controlled by a combination of both physical and geochemical parameters. Time scale analysis of O2 supply and consumption indicates that the occurrence of anaerobic hotspots is determined by the balance between characteristic times of O2 diffusion through the soil profile or exchange between macropores and the soil matrix, and the characteristic time of O2 consumption. Motivated by these findings, the present model was applied to simulate GHG cycling at a field site in Ontario, Canada, featuring macroporous agricultural soil. The model generally reproduced observed spatiotemporal variations in pore gas O2 and GHG concentrations and fluxes. Simulation results suggest that at this site, gas exchange processes between regions of preferential flow and the soil matrix play an important role in controlling the spatial variation of O2 in the soil. The simulations illustrate how the release of GHGs from the soil matrix into the macropores can lead to GHG emissions to the atmosphere. The modeling approach presented here, including dual domain flow and solute transport, as well as dual domain gas transport, shows promise for future studies with the aim to develop a more complete understanding of how soil structure affects complex GHG cycling in agricultural soils.

Groundwater with total dissolved sulphide concentrations in excess of 1.0×10−4 mol L−1 3 mg L−1 is relatively common at intermediate depths in sedimentary basins. However, the mechanisms responsible for the formation and spatial distribution of these sulphidic waters in sedimentary basins, which have been affected by periods of glaciation and deglaciation, are not fully understood. Sulphate reduction rates depend on many factors including redox conditions, salinity, temperature, and the presence and abundance of sulphate, organic matter, and sulphate-reducing bacteria. Two-dimensional reactive transport modelling was undertaken to provide potential explanations for the presence and distribution of sulphidic waters in sedimentary basins, partially constrained by field data from the Michigan Basin underlying Southern Ontario, Canada. Simulations were able to generally reproduce the observed depth-dependent distribution of sulphide. Sulphate reduction was most significant at intermediate depths due to anoxic conditions and elevated sulphate concentrations in the presence of organic matter in waters with relatively low salinity. The simulations indicate that glaciation-deglaciation periods increase mixing of waters at this interfacial zone, thereby enhancing rates of sulphate reduction and the formation of sulphide. In addition, the simulations indicate that glaciation-deglaciation cycles do not significantly affect sulphide concentrations in low permeability units, even at shallow depths (e.g., 25 m), while concentrations in permeable units remain stable below depths of 500 m.

Agricultural drainage ditches are necessary for regulating moisture contents in fields for crop production, and moreover, they can provide and regulate important ecosystem services. Drainage ditches and their associated riparian zones can be significant sources of enhanced N2O, CH4 and CO2 emissions, since they receive nutrient laden runoff and drainage from adjacent fields, while providing moisture conditions suitable for the production of greenhouse gases (GHGs). In this study, a uniform reactive transport model with up-scaled rate parameters was used to assess the dynamics of GHG cycling and biogeochemical processes in a crop field and riparian soils adjacent to a drainage ditch (i.e., located on the shoulder and midslope of the ditch) in eastern Ontario, Canada. Simulations adequately reproduced observed spatial variations in pore gas GHG concentrations in the cropped field and the monitoring locations near the drainage ditch. Spatial variations can be attributed to differences in environmental factors (i.e., soil moisture content and temperature), soil nutrient supply (i.e., NH4+, NO3–), and organic C availability. Simulation results also showed that at this site, GHG emissions did not vary significantly between monitoring locations, suggesting that the ditch slopes are not significant GHG emission hotspots. Both major rainfall events and warmer summer temperatures led to the development of hot moments for CO2 and N2O emissions and corresponding soil gas concentrations. Although the uniform model was successful in reproducing long-term trends in observed data, the model was not able to fully capture in-season temporal variability of fluxes and concentrations related to acute precipitation events (i.e., hot moments). Nevertheless, the modeling approach presented shows promise for future studies comparing the effect of management styles for agricultural soils and drainage ditch zones on GHG emissions and uptake.

The oxalate carbonate pathway (OCP) of the Iroko tree has been shown to function as a carbon sink, due to the precipitation of calcium (Ca) carbonate in the tree and the soil around its roots, a soil developed on a carbonate free bedrock. A 1D reactive transport (RT) model of the Iroko OCP system was constructed and used in a sensitivity analysis to replicate the geochemical pathway of the OCP and carbon sequestration, constrained by the Ca mass balance. The sensitivity analysis examined ranges of contributions for Ca sources to the system, including rainfall, dust, plant decomposition, mineral weathering, and optionally groundwater Ca sources. The study focused on identifying possible scenarios by which the observed magnitude of accumulated soil calcite could be reproduced, related to the OCP geochemical pathway. In line with field observations of the Iroko, all model results showed a soil transition from acidic (pH 5.1) to mildly alkaline (maximum pH 8.1) conditions, as well as precipitation of calcite. The location of calcite precipitation remained within the upper reach of the soil column, and between 29% and 63% of carbon was retained as calcite. The magnitude of calcite precipitation was limited both by the availability of Ca from the various sources, and by the tree's rate of root water and solute uptake. The sensitivity analysis confirms that surficial sources such as rainfall and dust deposition contribute only a small fraction of Ca, and indicate that mineral weathering and the extent of the root system, must play a dominant role in supplying Ca to facilitate carbon sequestration through calcium carbonate formation.

We investigate the effect of ice sheet geometry on groundwater flow patterns and meltwater ingress in a hypothetical sedimentary basin. The simulation results indicate that meltwater ingress is much greater in 3D domains with a relatively narrow ice sheet extent compared to a wide ice sheet, or a simplified 2D model. In high permeability units (HPUs), the simulated meltwater penetration depth can reach up to 750 m in 3D domains compared to 400 m in a comparable 2D domain. In low permeability units (LPUs), very limited meltwater penetration occurs, indicating that ice sheet geometry and model dimensionality do not substantially affect water flow in these units. A conservative tracer, with a source located at a depth of 500 m in both HPUs and LPUs, illustrates that solutes can be transported to greater depths in HPUs for a narrow ice lobe scenario, in comparison to a wide ice sheet or the 2D approach. Tracer transport in LPUs is unaffected by ice sheet geometry. Similarly, simulation results indicate that the ingress of dissolved oxygen (O 2 ) into HPUs is most substantial below a narrow ice lobe, while O 2 ingress into LPUs is not affected by ice sheet geometry. The numerical experiments indicate that 3D analysis will give more comprehensive results for flow patterns and reactive solute transport subjected to glaciation/deglaciation cycles in the case of a narrow ice lobe, but also suggest that a 2D approach might provide an adequate representation for the case of a relatively wide ice sheet.

Dipping anisotropy is a common feature in heterogeneous porous media that can substantially affect solute transport. For problems with complex geometry the influence of dipping anisotropy must be analyzed using numerical models, since suitable analytical solutions are not available. The most straightforward approach is to use a Cartesian coordinate system aligned with the material coordinate system. However, this approach is usually not practical, especially in 3D simulation domains with dipping layers and heterogeneous material properties. Furthermore, in the case of diffusion-dominated transport, the effect of anisotropy is often neglected. In this research, a general-purpose, fully 3-D unstructured grid code was developed to simulate diffusion-dominated solute transport in systems with dipping anisotropy, while accounting for complex geometry. The code has been verified against both 2-D and 3-D analytical solutions and has then been applied to two anisotropic diffusion problems, including an in-situ diffusion experiment and a hypothetical deep geologic repository, respectively. The simulation results indicate that consideration of anisotropy is required if the solute distribution in the rock matrix is of importance, in particular for assessing long-term evolution in layered systems. The formulation presented provides a versatile method for assessing diffusion-dominated solute transport in systems with dipping anisotropy subject to complex geometry.

The Canadian concept for a deep geological repository for used nuclear fuel includes highly-compacted and homogenized bentonite (HB) and low-alkali concrete such as Low-Heat High-Performance Concrete (LHHPC) as engineered barrier materials. The initial mineralogy and pore water compositions are different for bentonite, LHHPC and potential host rocks, such as granite and limestone. Consequently, chemical alterations are expected at the interfaces between these materials. Two scenarios were simulated to investigate the alterations at the material interfaces: 1) HB/LHHPC/host rock; 2) HB/host rock. For both scenarios, excavation damage zones (EDZs) were taken into consideration. Because of the low reactivity of bentonite compared to LHHPC, for the second scenario, only relatively minor mineral volume fraction and porosity changes were predicted for the HB/ granite and HB/limestone host rock cases, with limited impact on radionuclide migration. For the first scenario; however, due to the relatively high pH of pore water in LHHPC, substantial mineral dissolution and precipitation were predicted to occur at the HB/LHHPC and/or LHHPC/host rock interfaces, leading to porosity reduction or even pore clogging in close proximity (< 1 cm) of the interfaces. The predicted geochemical evolution depends mainly on the mineralogy and pore water chemical composition of the host rock. In the case of granitic host rock, Calcium Silicate Hydrate (C-S-H) phases initially present in LHHPC are predicted to transform into tobermorite, phillipsite, saponite and gypsum within 1,000 years. The simulations indicate a complex evolution of porosity, with an initial reduction, followed by a slight increase, and subsequent pore clogging. In the case of limestone host rock, saponite and sepiolite are the dominant minerals formed in the LHHPC. A sensitivity analysis shows that higher initial effective diffusion coefficients in the LHHPC enhance the diffusive fluxes of ions and result in earlier pore clogging in the host rock. The simulations illustrate that reactive transport modelling accounting for barrier reactivity and pore clogging provides a useful tool for assessing the migration of radionuclides across material interfaces for various canister failure scenarios.

In order to reduce contaminant mass loadings, thermal cover systems may be incorporated in the design of waste rock piles located in regions of continuous permafrost. In this study, reactive transport modeling was used to improve the understanding of coupled thermo-hydrological and chemical processes controlling the evolution of a covered waste rock pile located in Northern Canada. Material properties from previous field and laboratory tests were incorporated into the model to constrain the simulations. Good agreement between simulated and observational temperature data indicates that the model is capable of capturing the coupled thermo-hydrological processes occurring within the pile. Simulations were also useful for forecasting the pile's long-term evolution with an emphasis on water flow and heat transport mechanisms, but also including geochemical weathering processes and sulfate mass loadings as an indicator for the release of contaminated drainage. An uncertainty analysis was carried out to address different scenarios of the cover's performance as a function of the applied infiltration rate, accounting for the impacts of evaporation, runoff, and snow ablation. The model results indicate that the cover performance is insensitive to the magnitude of recharge rates, except for limited changes of the flow regime in the shallow active layer. The model was expanded by performing an additional sensitivity analysis to assess the role of cover thicknesses. The simulated results reveal that a cover design with an appropriate thickness can effectively minimize mass loadings in drainage by maintaining the active layer completely within the cover.

Placement methods and material availability during waste rock pile (WRP) construction may create significant heterogeneities in physical and geochemical parameters (such as grain size, permeability, mineralogy, and reactivity) and influence the internal pile structure. Due to the enormous scale of WRPs, it is difficult to capture the influence of heterogeneities on mine drainage composition and evolution. Although laboratory- or field-scale experimental studies have provided much insight, it is often challenging to translate these results to full scale WRPs. This study uses a numerical modeling approach to investigate the influence of physical and chemical heterogeneities, structure, and scale on the release of acid rock drainage (ARD) through 2D reactive transport simulations. Specifically, the sensitivity of drainage quality to parameters including grain size distribution, sulfide mineral weathering rates, abundance and distribution of primary minerals, and pile structure as a function of construction methods are investigated. The geochemical model includes sulfide oxidation, pH buffering by calcite dissolution, and ferrihydrite and gypsum as secondary phases. Simulation results indicate that the implications of heterogeneity and construction method are scale-dependent; when grain size distribution trends observed in a pile's core are applied to the entirety of a pile, results between push- and end-dumping methods vary substantially—however, predicted drainage for different construction methods become more similar when features such as traffic surfaces, structural variation, and multiple benches are also considered. For all scales and construction methods investigated, simulated results demonstrate that pile heterogeneity and structure decrease peak mass loading rates 2 to 3-fold, but cause prolonged ARD release compared to the homogeneous case. These findings have implications for the economics of planning water treatment facilities for life of mine and closure operations.

The oxidation of sulfide minerals such as pyrite present in waste rock results in elevated sulfate, enhanced metal loadings and in many cases low pH conditions. Recently, many mines have opened in remote areas, including regions subject to permafrost conditions. In these regions, freeze-thaw cycles and the possible development of permafrost in mine waste add to the complexity of weathering processes, drainage volumes and mass loadings. To assess weathering in these waste rock piles, the reactive transport code MIN3P-HPC has been enhanced by implementing constitutive relationships related to freeze-thaw cycles that control flow patterns, solute transport, generation and transport of heat, as well as geochemical reactions and their rates. Simulations of a hypothetical pyrite-rich waste rock pile placed onto natural permafrost were conducted under reference climate conditions. Additionally, the effect of a warming climate was also studied through a sensitivity analysis. The simulation results indicate a potentially strong coupled effect of sulfide mineral weathering rates and a warming climate on the evolution and persistence of permafrost within waste rock piles and the release of acidic drainage. For relatively low sulfide mineral oxidation rates, the simulations indicate that permafrost can develop within waste rock piles, even under warming climate conditions. However, the results for low reactivity also show that mass loadings can increase by >50% in response to a slight warming of climate (3°C), relative to reference climate conditions. For the chosen reference reaction rates, permafrost develops under reference climate conditions in the simulated waste rock pile; however, permafrost cannot be maintained for a marginally warmer climate, leading to internal heating of the pile and substantially increased production of acidic drainage (>550%). For high reaction rates, the simulations suggest that internal heating takes place irrespective of climate conditions. Evaluation of thermal covers indicates that significant reductions of mass loadings can be achieved for piles with low and reference reactivity (91–99% in comparison to uncovered piles), but also suggest that thermal covers can be ineffective for piles with high sulfide content and reactivity. Together, these simulations provide insights into the complex interactions controlling waste rock weathering in cold-region climates.

The numerical simulation of flow and reactive transport in porous media with complex domains is nontrivial. This paper presents a method to implement fully unstructured grid capabilities into the well-established software ParMIN3P-THCm, a process-based numerical model designed for the investigation of subsurface fluid flow and multicomponent reactive transport in variably saturated porous media with parallelization capability. The enhanced code, MIN3P-HPC, is modularized to support different cell types, spatial discretization methods and gradient reconstruction methods. MIN3P-HPC uses a vertex-centered control volume method with consideration of both vertex-based and cell-based material properties (e.g., permeability). A flexible parallelization scheme based on domain decomposition and thread acceleration was implemented, which allows the use of OpenMP, MPI and hybrid MPI-OpenMP, making optimized use of computer resources ranging from desktop PCs to distributed memory supercomputers. The code was verified by comparing the results obtained with the unstructured grid version to those produced by the structured grid version. Numerical accuracy was also verified against analytical solutions for 2D and 3D solute transport, and by comparison with third-party software using different cell types. Parallel efficiency of OpenMP, MPI and hybrid MPI-OpenMP versions was examined through a series of solute transport and reactive transport test cases. The results demonstrate the versatility and enhanced performance of MIN3P-HPC for reactive transport simulation.