Accurate simulation of CO2 sequestration in deep saline aquifers requires consistent treatment of history-dependent processes. Operator-Based Linearization (OBL) provides an efficient and robust framework for compositional flow simulation. However, it was originally designed for reversible physical processes, since the operator evaluation depends solely on the instantaneous thermodynamic state. In this work, we extend the OBL approach by introducing history-dependent variables into the operator parameter space while keeping them outside the Newton unknown vector, thereby augmenting OBL with history dependency while preserving the dimension of the Jacobian matrix. To apply the augmented OBL framework to hysteresis modeling, a new hysteresis algorithm is proposed, in which the maximum gas saturation is used as the history coordinate; the flow equations are solved fully implicitly with respect to the primary variables, while this history variable is held fixed during Newton iterations and updated locally after timestep convergence to capture drainage–imbibition transitions and CO2 dissolution feedback. The augmented OBL approach is validated against the academic DARSim simulator and the commercial CMG simulator, demonstrating close agreement with conventional hysteresis modeling approaches. A sensitivity analysis of the operator parameterization shows that the resolution of the parameter space governs both accuracy and computational cost within the OBL framework. Numerical experiments in one-dimensional homogeneous and two-dimensional heterogeneous models further demonstrate that the proposed approach provides a practical and accurate means of incorporating history-dependent physics into reservoir-scale compositional simulations of CO2 sequestration.

The effective management of geoenergy systems heavily relies on robust modeling frameworks that integrate diverse simulation capabilities, including flow and transport, phase equilibrium, geochemistry, and geomechanics. While a multiphysics simulation engine within a unified framework has its advantages, integrating specialized modeling packages often enhances viability. Efficient and seamless communication between these engines becomes crucial for improving the performance and scalability of the integration. Advanced parametrization techniques can facilitate this integration by efficiently approximating and interpolating coupling data, ensuring both speed and accuracy. In this study, we compare the efficiency of different interpolation techniques used for the parametrization of complex many-component fluid systems in compositional simulation. We use an operator-based linearization (OBL) framework that leverages the general formulation of the corresponding conservation laws. OBL effectively learns the operators required for the assembly of these laws, while interpolation delivers fast evaluation of operators and their derivatives for all physical states in a simulation domain. Multilinear interpolation is a simple and robust approach, yet it has poor scaling properties with respect to the dimension of the physical state. To alleviate interpolation costs in multiple dimensions, we study the performance and accuracy of other interpolation techniques, including linear interpolation with standard and Delaunay triangulation. Overall, this approach provides great flexibility, saves development costs, and simplifies the incorporation of thermodynamics and geochemistry engines for precise modeling of phase equilibrium, reactive transport, dissolution-precipitation, and kinetics of chemical reactions. We investigate the efficiency of the aforementioned interpolation strategies in terms of robustness, accuracy, performance, and scalability. As test models, we consider a multicomponent fluid model with thermal and reactive effects relevant to carbon dioxide (CO2) sequestration applications. As expected, we observe the performance gain of linear interpolation compared with multilinear interpolation, increasing with the number of components. We investigate the robustness and efficiency of standard and Delaunay triangulations for systems with a significant number of components. This research extends the scalability of the OBL framework and addresses the challenges of high dimensionality in compositional modeling. Consequently, this approach holds significant potential for integrating various complex multiphysics problems, advancing the development of more comprehensive digital twins for geoenergy systems management.

Injection of \ce{CO2}-saturated water into carbonate rocks drives reactive transport in which mineral dissolution alters porosity and permeability, focuses flow, and can organize into wormholing dissolution patterns. Predicting these patterns and the associated breakthrough behaviour requires a modelling framework that resolves the full coupling between multiphase flow, multicomponent transport, and geochemistry across dissolution regimes, and that remains computationally feasible at field-relevant resolutions. Building on the operator-based linearization (OBL) framework for fully implicit coupling of compositional reactive transport and external geochemical solvers developed in the companion paper \citep{novikov2026partI}, we apply the coupled formulation to heterogeneous two- and three-dimensional simulations of carbonated-water injection into limestone. Geochemical equilibrium and speciation are evaluated by PHREEQC and Reaktoro and incorporated as state-dependent operators, while the global flow and transport problem is solved fully implicitly within open-DARTS or GEOS. The simulations reveal a pronounced grid-resolution dependence of the computed dissolution structures: finer grids produce thinner, more sharply defined wormholes and earlier breakthrough, a trend confirmed by ensemble studies over multiple heterogeneous realizations and by fully three-dimensional core-scale simulations on unstructured meshes. The porosity--permeability relationship strongly controls channel morphology and breakthrough efficiency, while the dependence of pore volumes to breakthrough on injection rate is weakly non-monotonic and remains within a wormhole-dominated regime. A scalability study with up to millions of cells shows that the OBL treatment of geochemical coupling is not the dominant computational cost; the principal expense shifts to the global linear solver and preconditioner at high parallel counts. The results demonstrate that the proposed framework provides a scalable basis for high-fidelity simulation of reactive dissolution, while highlighting the intrinsic multiscale nature of wormholing and the need for grid-independent modelling strategies.

Thermal-Hydro-Mechanical-Compositional analysis is crucial for addressing challenges like wellbore stability, land subsidence, and induced seismicity in the geo-energy applications. Numerical simulations of coupled thermo-poromechanical processes provide a general-purpose tool for evaluating these phenomena across laboratory and field scales. However, efficient integration of the coupled equations for fluid mass, energy and momentum poses multiple numerical and implementation difficulties, such as combining different numerical methods on staggered grids and associated limitations on admissible grids. This paper introduces a novel fully-implicit Finite Volume Method (FVM) for modeling thermal compositional flow in thermo-poroelastic rocks. The scheme employs gradient-based, coupled multi-point approximations of fluid mass, momentum and heat fluxes.The novelty of the scheme lies in its integration of temperature as a parameter in the flux approximation process. The scheme supports a wide range of cell topologies, arbitrary heterogeneity and anisotropy as well as various boundary conditions, while respecting local flux balance under temperature gradients. Overall, the scheme represents a unified FVM-based approach for the integration of all conservation laws relevant to geo-energy applications on a cell-centered collocated grid. Additionally, the implemented two-stage block-partitioned preconditioning strategy enables the efficient solution of obtained linear systems.The framework, implemented in the open-source Delft Advanced Research Terra Simulator (open-DARTS), leverages the Operator-Based Linearization (OBL) technique for flexibility in compositional fluid properties. Rigorous validation demonstrates the framework’s capabilities in capturing advanced phenomena, including thermal expansion, thermo-poroelastic effect and compositional flow with phase transitions. The performance of preconditioning strategy is assessed using the mechanical extension of the SPE10 benchmark model.

Abstract Pore pressure fluctuation in subsurface reservoirs and its resulting mechanical response can cause fault reactivation. Numerical simulation of such induced seismicity is important to develop reliable seismic hazard and risk assessments. However, modeling of fault reactivation is quite challenging, especially in the case of displaced faults, i.e., faults with non-zero offset. In this paper, we perform a systematic benchmarking study to validate two recently developed numerical methods for fault slip simulation. Reference solutions are based on a semi-analytical approach that makes use of inclusion theory and Cauchy-type singular integral equations. The two numerical methods both use finite volume discretizations, but they employ different approaches to represent faults. One of them uses a conformal discrete fault model (DFM) while the other employs an embedded (non-conformal) fault model. The semi-analytical test cases cover a vertical frictionless fault, and inclined displaced faults with constant friction and slip-weakening friction. It was found that both numerical methods accurately represent pre-slip stress fields caused by pore pressure changes. Moreover, they also successfully cope with a vertical frictionless fault. However, for the case with an inclined displaced fault with a constant friction coefficient, the embedded method can not converge for the post-slip phase, whereas the DFM successfully coped with both constant and slip-weakening friction coefficients. In its current implementation, the DFM is therefore the model of choice when accurate simulation of local faulted systems is required.

Simulator (Voskov et al., 2023) is a simulation framework for forward and inverse modelling and uncertainty quantification of multi-physics processes in geo-engineering applications such as geothermal, CO2 sequestration, water pumping, and hydrogen storage.To efficiently achieve high levels of accuracy on complex geometries, it utilizes advanced numerical methods such as fully implicit thermo-hydro-mechanical-chemical (THMC) formulation, a highly flexible finite-volume spatial approximation and operator-based linearization for nonlinear terms.open-DARTS goals are computational efficiency, extensibility, and simplicity of use.For this reason, open-DARTS is based on a hybrid design with an efficient core C++ implementation wrapped around a highly customizable and easy-to-use Python code.

Direct numerical simulations of the evolution of turbulent spots in the boundary layer on a flat plate at a zero angle of attack at the freestream Mach number M ∞ = 6 are carried out. The propagation of artificially excited localized three-dimensional vortical disturbances with different initial amplitudes, which, when propagating downstream, develop into turbulent spots, is considered. A direct numerical simulation is performed by solving the Navier–Stokes equations for three-dimensional compressible gas flows using the in-house solver that implements an implicit shock capturing numerical scheme. It is shown that the universal quasi-monotonic numerical scheme makes it possible to correctly estimate the main characteristics of turbulent spots: the transverse spreading angle and the velocities of the leading and trailing fronts. The agreement between the parameters of the obtained spots and the results of other authors is demonstrated.

Elliptic differential operators describe a wide range of processes in mechanics relevant to geo‐energy applications. Extensively used in reservoir modeling, the Finite Volume Method with TPFA can be consistently applied to discretize only a specific type of application under severe assumptions. In this paper, we introduce a positivity preserving Nonlinear Two Point Stress Approximation (NTPSA) based on the recently developed collocated Finite Volume scheme for linear elastic mechanics. The gradient reconstruction is different from the one used in Nonlinear TPFA, but a similar form of weighting scheme is employed to reconstruct the traction vector at each interface. The convergence of the scheme is tested with a homogeneous anisotropic stiffness tensor. The motivation behind the implementation of a new discretization framework in mechanics is to develop a uniform discretization technique preserving monotonicity for generic poromechanics applications.

ABSTRACT Quantification of the poromechanical response of subsurface formations due to human-induced pore pressure fluctuations is critical for the performance and stability assessment of many geo-energy systems. In particular, natural faults in the subsurface introduce the hazard of induced seismicity. Numerical modeling of fault reactivation is challenging, while the specific details of induced stresses and fault slip in reservoirs with displaced (i.e. non-zero offset) faults may cause additional challenges depending on the type of numerical formulation employed. To facilitate the systematic development and testing of numerical tools for the simulation of induced seismicity in faulted reservoirs we developed a set of semi-analytical test problems of increasing complexity, based on inclusion theory and Cauchy singular integral equations. With these we investigate the accuracy of two recently developed Finite Volume (FV) schemes with collocated and staggered arrangements of unknowns. One of them employs a conformal discrete fault model (DFM) which can guarantee sufficient accuracy at the cost of adaptive mesh refinement but may suffer from modelling and computational challenges when addressing large-scale realistic geological configurations. The second one employs an embedded (or non-conformal) discrete fault model (EDFM) which avoids the need for excessive mesh refinement, but of which the accuracy and the range of applicability are still to be investigated. We found that both numerical schemes accurately represent the pre-slip Coulomb stresses, but show different degrees of accuracy in representing the resulting depletion-induced fault slip. The semi-analytical benchmark data are available via DOI 10.4121/22240309. INTRODUCTION The kernel of this paper is formed by a series of semi-analytical poro-mechanical test problems of increasing complexity with the aim to systematically compare the capacities of two poro-mechanical finite-volume-based simulation codes: one developed by Novikov et al. (2022b) which employs a discrete fault model (DFM). It forms part of a comprehensive porous media simulation package, the Delft Advanced Reservoir Terra Simulator (DARTS) and will be referred to with that acronym. The second code, developed by Shokrollahzadeh Behbahani et al. (2022), is based on a smoothed version of the embedded discrete fault model (sEFVM) and will be referred to with that last acronym. Both codes are being developed as part of the DeepNL Science4Steer project (NWO, 2017), and Appendices A and B give a brief overview of their characteristic features. The grids used in this study are presented in Appendix C.

The problem of acid treatment in carbonate reservoirs at different scales is considered. The problem of two-phase multicomponent Darcy flow through a porous medium with chemical reactions is solved at the core scale. The regimes of acid injection into the core that lead to the wormholing phenomena are investigated. The present study makes it possible to determine the wormhole growth rate at various reaction kinetics parameters and injection rates and for various permeability distributions. The dependences obtained are used in the large-scale model of fracture acidizing in a reservoir (FAR) that also includes the model of two-phase multicomponent flow with chemical reactions through a porous medium and the model of acid transport through a fracture. Additional contribution of wormholes to the conductivity is simulated. Based on the model proposed, a series of calculations of fracture acidizing in a reservoir is carried out, the effect of taking into account the wormholes in the large-scale model proposed is investigated, and the influence of the constitutive parameters of the process considered such as the fracture length, the reservoir permeability, and the solution injection rate on the effectiveness of treatment estimated in terms of the production ability of well after treatment is studied.

We present a scalable collocated Finite Volume Method (FVM) to simulate induced seismicity as a result of pore pressure changes. A discrete system is obtained based on a fully-implicit fully-coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a flexible unstructured mesh. The cell-centered collocated scheme leads to a convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a generic multi-point flux approximation is formulated to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable processing algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate and efficient; and provides a promising framework for accurate and efficient simulation of induced seismicity in various geoscientific applications.

Changes in temperature and pressure conditions during the development of an oil reservoir with a high paraffin content in a limiting saturated state causes an in-situ phase transition to a solid state, and the filtering of a mixture of oil with solid paraffin particles at a temperature below the flocculation temperature causes clogging in the pore space of the reservoir in narrow places and bottle necks of pores. Thus, laboratory studies are carried out to determine the critical points and parameters of an oil - gas - paraffin mixture under various thermobaric conditions. A mathematical model is developed that describes wax deposition in the pore space of a low-temperature oil reservoir and makes it possible to calculate its permeability during development. The model parameters are adapted to their experimental values. Numerical calculations are used to determine the values of the main parameters that affect the reservoir permeability.

In this work, numerical studies are conducted for a supersonic M = 2.0 steady flow over a flat plate with a suction through single circular hole of variating intensity. The computed flow fields near the suction hole are estimated to build special boundary condition of Dirichlet type. Such a simulated suction can be used in DNS of disturbances propagation over walls with multiple suction holes.

Numerical and theoretical stability studies are conducted for a supersonic M = 2.0 flow over a flat plate with a spanwise strip of distributed suction. Propagation of artificially excited 3D wave trains relevant to the first-mode instability are considered. Direct numerical simulations (DNS) as well as assessments based on the linear stability theory (LST) are performed using the in-house code “HSFlow++” (High Speed Flow solver). It is shown that, despite strong spatial nonuniformities of the mean flow near the suction strip ends, the e-N method based on the local-parallel LST gives reasonable estimates of the laminar flow control performance of the suction strip.

Low-thrust rocket engines are widely used in rocket and space technology for correcting the position of a spacecraft in orbit, for controlling motion along a trajectory, etc. Their number in the propulsion system can be from one to tens of units. Accordingly, the efficiency of their work significantly affects the perfection of the propulsion system as a whole. The object of the study was the low-thrust rocket engine combustion chamber operating according to the gas-liquid scheme. There were performed computational and parametric studies of various factor effects on the characteristics of the working process in the combustion chamber. The dependences of the coefficient of the consumable complex and parameters of the working process of the low-thrust rocket engine chamber on the influencing factors when using ethanol and kerosene as a fuel were calculated. A comparative analysis of the results of using these two components under similar conditions was carried out, which made it possible to reveal the influence of the physicochemical properties of the combustible component on the efficiency of the working process organization. The results obtained can be used in the design of low-thrust engines operating on the kerosene–oxygen and ethanol–oxygen propellants.

We develop a collocated Finite Volume Method (FVM) to study induced seismicity as a result of pore pressure fluctuations. A discrete system is obtained based on a fully-implicit coupled description of flow, elastic deformation, and contact mechanics at fault surfaces on a fully unstructured mesh. The cell-centered collocated scheme leads to convenient integration of the different physical equations, as the unknowns share the same discrete locations on the mesh. Additionally, a multi-point flux approximation is formulated in a general procedure to treat heterogeneity, anisotropy, and cross-derivative terms for both flow and mechanics equations. The resulting system, though flexible and accurate, can lead to excessive computational costs for field-relevant applications. To resolve this limitation, a scalable parallel solution algorithm is developed and presented. Several proof-of-concept numerical tests, including benchmark studies with analytical solutions, are investigated. It is found that the presented method is indeed accurate, stable and efficient; and as such promising for accurate and efficient simulation of induced seismicity.

The authors have analyzed the essence, structure, rates of negative trends in the Arctic coastal zone, risk factors, the consequences of climate change. They have outlined ways to mitigate and overcome the trends, paying particular attention to the elimination of objects of accumulated environmental damage in the Russian Arctic zone. The article clarifies the concept and methods of assessing climate risk. The researchers have carried out a PEST-analysis for the factors of influence of the external environment on the infrastructure formation of the Arctic coastal territories and proposed to use the mechanism of public-private partnership, develop renewable energy sources, conduct ethnological expertise of projects, create special units for managing climate risks in stakeholder companies.

The article studies issues of land tenure planning for implementation of projects aimed at industrial development of the Arctic. Using the example of Northern provinces of Canada it shows evolution of land tenure strategic planning, analyzes its role in social and economic development of the territory. It is shown that involvement of aboriginal people of the North in the process of planning the use of land, forest and other natural resources can lower conflicts among land users, mining companies and the local population, protect territories of traditional land tenure in places of residence and traditional natural resource use of aborigine people and create necessary conditions for the development of traditional types of activity and sustainable space development of the Arctic. Canadian experience of land tenure planning in development of Arctic territories in the area of aboriginal people residence can be used in the Arctic zone of the Russian Federation to balance interests of concerned parties, i.e. local bodies of power, business and aboriginal people of the North.

Acid impact is the highly efficient kind of carbonate reservoir stimulation. In contrast to hydraulic fracturing acid impact improves matrix permeability that remains high while production. In this case propagation of reactive flow into matrix brings stimulation that means that the flow factor is crucial for fracture acidizing. This paper is focused on the reactive flow modelling in the existing fracture and oil saturated matrix. Two phase multicomponent Darcy flow which takes into account dissolution kinetics is modelled in two domains (fracture, matrix). Using prescribed porosity permeability relationship the stimulated permeability field is calculated. Acid propagation into matrix is investigated for different injection control regimes. Stimulated productivity is assessed at the final stage of stimulation in the enlarged domain. Obtained results show that the injection time plays a key role for stimulated productivity. Sensitivity analysis to the other parameters is also given.