Reactive transport in porous media: a review of recent mathematical efforts in modeling geochemical reactions in petroleum subsurface reservoirs

The rapid advancements in the computational abilities of numerical simulations have attracted researchers to work on the area of reactive transport in porous media to improve the hydrocarbon production processes from mature reservoirs. In the hydrology community, reactive transport is well developed where the main research focuses on studying the movement of groundwater and contaminants in aquifers, and quantifying the effect of chemical reactions between the rocks and water. Recently, great efforts have been made to adapt similar models for petroleum applications where multiphase flow is experienced in the subsurface reservoirs. In such systems, thermodynamic and chemical equilibrium is key in establishing an accurate description of the states of the fluids existing in the reservoir. This paper presents a detailed and comprehensive review on the concepts of geochemical modeling, and how it can be mathematically adapted to petroleum multiphase flow problems in porous media. We introduce key physical concepts outlining the treatment of chemical reactions in experimental trials and then explain in detail how a network of chemical reactions can be modeled mathematically for numerical simulation applications. The steps of characterizing the physical behavior of the fluid flow—through a set of governing equations by either natural or molar variables formulations, and the methodology to simplify and incorporate the numerical algorithms into an existing reservoir simulation scheme are shown as well. We finally present two numerical cases from the literature to highlight the key variations between the different variable formulations and comment on the advantages and disadvantages of each approach.


Introduction
The topic of reactive transport has attracted the efforts of researchers from different disciplines for the past four decades [17,21,26,27,73,98,99] due to its role in understanding the composition of natural water aquifers and the formation and dissolution of rocks and minerals in geologic formations such as deep hydrocarbon reservoirs. The general framework of a reactive transport model was mainly developed to understand the rock-fluid interactions in groundwater movement. At early stages, the modeling process was restricted to single components in the hydrology community, but was later extended to include multicomponents after Litchner established the theoretical basis of the continuum model [73]. This introduced new capabilities to treat the coupling of various transport and biochemical processes. Such an advancement paved the path for many alterations of the original framework to be adapted for reactive contaminant transport, non-isothermal water-rock interactions with kinetic reactions and many more [14,25,77,100,120].
Traditionally, the solution schemes in numerical codes for reactive transport problems can be divided into two families [90]. The first consists of codes that utilize the sequential approach and can be broken down into two main subcategories: sequential iterative (SIA) and sequential non-iterative (SNIA). In the SNIA approach, the mass transport and chemical reactions are decoupled and solved separately where constant parameters are used when computing the pressure gradient. This is referred to as operator splitting, where the flow equations are solved at each time step without the reaction terms and then followed by the reaction step. This approach introduces splitting errors that is usually reduced by utilizing the SIA. The latter approach is quite similar in the fact that the mass and chemistry equations are decoupled, but the mole fractions and phase velocities from the transport step are used in a chemical submodel to update the concentrations of the components while iterating between solutions of transport and chemistry until an error tolerance is met. The operator splitting techniques have been utilized in many simulators such as HYDROGEOCHEM [122], THOUGHREACT [121], PARSIM [15] and PHREEQC [28], as these approaches are characterized by their easy implementation in a working simulator and fewer demands of computational resources in terms of memory and CPU time. Moreover, legacy codes can be coupled for each phenomenon individually since both chemistry and transport are treated as operators that are modeled independently of one another. However, the SNIA and SIA require small time steps to be used, which could lead to slow convergence rates. This issue can be overcome by introducing the global implicit approach (GIA), which we classify as the second family of solution techniques and is mainly utilized by CrunchFlow [19], MIN3P [101] and PFLOTRAN [56] simulators. For simple models with a limited number of chemical components, the fully implicit/coupled integration method allows solving the system directly. The scheme is often accurate as both the transport and chemical steps are solved simultaneously using a Newton-Raphson numerical approach with less restrictions on time step size, providing stability and robustness [42]. However, if the number of components in the system is large, this method suffers from an excessive computation time and memory requirements, especially in twoand three-dimensional simulations. The difference in the time scales of the transport and the chemistry steps and degree of coupling might impose additional hurdles to the efficiency of the approach. A summary of the properties of the GIA and SIA approaches is presented in the flowchart in Fig. 1, and a detailed comparison of key numerical capabilities of reactive transport modeling codes can be found in the work of [100].
In groundwater simulations, many codes pioneer in modeling reactive transport on the basis of the continuum model, such as HYDROGEOCHEM [122], TOUGHREACT [121], MIN3P [101], PFLOTRAN [56] and CrunchFlow [19]. A famous package is the one developed by Lawrence Berkeley National Laboratory, TOUGHREACT [121], and has been tested extensively for problems with multiple phases and various rock-fluid interactions. It mainly employs the sequential iteration scheme discussed earlier and can handle thermal conditions in the reservoir. Another package is CrunchFlow [19] that utilizes fully implicit schemes for 1-D and 2-D problems but lacks the treatment of multiphase flow. However, the PFLOTRAN package developed by Litchner [56] can deal with multicomponent and multiphase systems but is limited to a defined set of problems. The phase behavior is fully characterized using Henry's law approach, and it supports parallel computations. Although these sets of codes are relatively mature, they cannot be applied directly to oil production scenarios because of the absence of the hydrocarbon phase treatment in the system.
On the other hand, and in the context of petroleum reservoir simulation, very limited models have been developed to couple transport and chemistry for hydrocarbon systems. Amongst these packages are STARS [53], ECLIPSE 300 [76] and UTCHEM [83], where the first two are commercialized, while the third serves as an open-source research simulator. STARS [53] is a module incorporated in the Computer Modelling Group simulator, GEM [52], and it is a generalized thermal-compositional simulator with a full equation of state characterization for phase behavior coupled with mineral dissolution and precipitation kinetics. Similar features are observed in Schlumberger's simulator ECLIPSE 300 [76] with its robust performance in modeling the compositional behavior of multiphase systems. Both codes employ a fully implicit scheme in their numerical algorithms but lack comprehensive chemical reactions implementations. On the contrary, UTCHEM [83], a simulator developed at The University of Texas at Austin, offers full reactions characterization with abilities to model chemical flooding scenarios in hydrocarbon reservoirs with full thermodynamic capabilities. Despite these features, the simulator lacks in the area that only SIA is employed, which might cause accuracy problems to arise, especially for multicomponent and multiphase systems in petroleum reservoirs. A fully general GIA implementation of coupled transport and chemistry for petroleum reservoir simulation applications does not yet exist.
It is important to distinguish between reactive transport modeling and geochemical modeling as most people refer to both terms interchangeably. In geochemical modeling, the interactions between water and earth material are analyzed and interpreted in the light of the chemical reactions in the system to provide useful insights on the composition and the thermodynamic behavior of such studied systems. The analysis is usually carried out in a 0-D environment to focus on the effects of temperature, ion speciation, sorption and solubility on the concentration of the various minerals in the system. On the other hand, reactive transport modeling focuses on more complex physical processes involving fluid flow and mass transfer, which require establishing geochemical equilibrium and understanding the thermodynamical state of the system [57,58]. This major distinction challenged the abilities of traditional reactive transport models to deal with multiphase flow of fluids in heterogeneous porous media, especially in hydrocarbon reservoirs. When multiple fluids are present in the same space, determining the geochemical equilibrium of a multicomponent system and the phase split ratios can be gruesome. Thus, incorporating geochemical reactions with the transport of fluids in a hydrocarbon reservoir is not a straightforward process but rather an extensive procedure at the physical and numerical level.
Recent advancements were able to modernize the approach of reactive transport modeling, thus coupling the flow with a comprehensive thermodynamic Review Paper SN Applied Sciences (2021) 3:401 | https://doi.org/10.1007/s42452-021-04396-9 description of the system and the kinetic framework [43-45, 62, 75, 95, 108]. However, these models cannot be treated as a traditional black-box simulator. In fact, special treatments should always be applied to the mathematical formulation of the problem based on the physical application at hand. The skill of selectively filtering the processes that affect your given system and neglecting the impact of minor processes requires a deep understanding of the geological and chemical settings of the system. This essential understanding highlights the complexity of reactive transport modeling as the input of the system is highly dependent on sources and information from various fields including but not limited to chemistry, thermodynamics, fluid dynamics and rock properties. This in turn requires researchers and scientists to have high expertise in multiple fields to be able to fully characterize and describe the problem at hand. The objective of this work is to provide a central source describing previous key breakthroughs in the area of reactive transport modeling with a special focus on its applications in multiphase flow in porous media. We will give great attention to the different mathematical formulations used to couple flow and reactions in a reservoir simulation framework and critically analyze the performance of each approach. We will also present the current challenges facing the development of comprehensive geochemical reservoir simulations for different practical applications and finish our review with concluding remarks on the potential of the implementation of reactive transport simulators in the field.

Experimental reactive transport
Typically, a reactive transport system is studied quantitatively through modeling tools that rely on mathematical formulations of the problem. However, it is often challenging to account for complex reaction systems in a simple yet physically comprehensive mathematical formulation. This complexity in the modeling highlights the need to establish an experimental framework for controlled laboratory experiments to better identify the processes and reaction paths in a reactive transport system. The experimental results contribute massively in improving existing models through providing high-resolution data sets.
The main key in designing laboratory experiments for reactive transport is to ensure a deep knowledge and understanding of the conceptual model of the physical process. This involves studying the degree of heterogeneity in the porous media and how would state variables in the governing equations vary with the progression of the reactions [24,88]. A prime focus in this case would be to understand the physics of the flow at both the pore scale and the continuum (Darcy) scale.
Flow tests at laboratory scales with varying configurations are usually carried out to examine the interactions between rock-fluid systems. Core flooding tests are performed to acquire petrophysical data on the properties of the porous media through specific flow cells, pumps and sensors, see Fig. 2 for a core flooding setup by [60]. The fluids in this setup are injected into the sample core, while confining stress is applied. The difference in the pressure is recorded in the chamber, and back pressure is applied to regulate the pressure at the core outlet. In some cases, using actual rock samples in reactive transport experiments might yield inconsistent results due to the extreme heterogeneity present in the rock. Hence, analog materials with consistent properties are used to facilitate the fundamental understating of the physical processes under study as the fabricated gypsum specimens by [34] and the textured glass (non-reactive) with a polished, flat, transparent rock crystal by [37].
In some cases of core flooding when a back pressure regulator is used, it is possible to carry out a chemical analysis on the effluent stream. The experiments in [94] focused on the dissolution of Marly Dolostone in CO 2 rich brine where the concentration of the calcium ions, magnesium ions and inorganic carbon was measured through the effluent sampling from a back pressure regulator. It is also possible to monitor the effluent chemistry directly through sampling probes across the flow cells, and such experiments are generally performed in 1-D and 2-D domains. Although 3-D experiments represent the boundaries and the conditions of the reservoir to an accurate degree, they are rarely performed due to the associated operational complexity. The experimental design for continuous sampling of the fluid is composed of a flow chamber that is filled with a material to resemble the porous media under investigation. The chamber has a typical maximum length or width of 10 cm and is connected with a feed line that matches the inlet components of your system. The final mixture after the reaction process ceases is collected at the other side of the chamber. The sampling ports are needed to take samples to measure the distribution of the solute concentration spatially and temporally. The measurements are further used to cross-check and validate the accuracy of the numerical models. It is important to note that only a limited number of samples are viable to avoid disturbing the dynamics of the fluid flow in the cell, see Fig. 3.
The work of [69] presented a comprehensive investigation for diffusive and advective processes to observe the effect of porosity changes on the transport parameters. This was accomplished through injecting a reactive solution (acidic zinc sulfate) into the pore space causing gypsum and smithsonite formation. The experiments were later used to benchmark the reactive transport code Hytec [105]. Later came the work of [103] where he combined laboratory experiments with a deep investigation on the changes of rock properties at the pore level. The information obtained allowed for a better representation of the continuum-scale model and introduced specific mixing parameters to the advection-dispersion reaction equations to account for highly non-uniform pore-scale concentration gradients. Similarly, [64] performed a laboratory experiment in a 2-D flow cell with multiple sampling ports and demonstrated the shortcomings of the reactive continuum model to reproduce the concentration of the solute measured at the sampling ports.
However, selective sampling does not efficiently contribute to the understanding of how the reaction rate progresses with space and time [29,35,39,41,51,54,93]. Only specific systems where the reaction might cause the products to change colors might be measured, but this only provides a local understating of the rate and does not correlate with the change in properties of the porous medium. To counter such a shortcoming, a computer-controlled triaxial system (CCTS) was designed and built in  MIT Geotechnical Engineering Laboratory which provides a continuous effluent chemistry monitoring system. The initial CCTS development of the system was inspired by the setup of [68] in which he designed a triaxial pedestal with integrated sensors for effluent chemistry analysis. The CCTS setup was progressed and updated over the years by [10,72,102], and it is extremely efficient in providing long-term continuous measurements of temperature and concentration and on the back-pressurized effluent stream, see Fig. 4.
In brief, laboratory experiments are essential to gain a deep understating of the conceptual and numerical model of the reactive transport process. However, the experimental setup can be challenging in many cases and upscaling problems might arise when mapping the data from flow cells to field-scale systems. Moreover, geochemical reactions are mainly controlled by processes occurring at the pore scale, which might alter the connectivity of the porous media and substantially affect the flow of fluids. The evolution of the porous media especially in heterogeneous rocks poses a challenge in representing reactive transport in a continuum approach. Traditionally, reactions and transport have been described on the basis of empirical relationships [104] and thus severely limiting the simulation results due to the loss of information in the upscaled parameters. Recent efforts have been utilizing the technological advancement in the area of CT-scanning and X-ray micro-tomography to provide a greater focus on the fluctuations of such parameters at the pore level of the rock.
The work of [18,38,82,107] has succeeded in improving pore-scale modeling formulation thankfully due to the increasing computational power of the computer machines, allowing for processing of very accurate and detailed images of the pore network that are obtained experimentally [7,8,89]. Although multiple efforts have been made in the recent years to advance pore-scale experiments, the process of upscaling the findings from a pore-scale level to a continuum level that extends in order of km's is still not completely feasible [59,84,87,115]. This is mainly attributed to the complicated experimental procedure and low reproducibility.

Assessing the role of CO 2 in geochemical reactions
Typically, crude oil production goes through three stages. The first stage is the primary recovery where the oil is ousted into the wellbore from the reservoir using its existing energy, such as compaction, gravity drainage, gas drive and water drive mechanisms. The next recovery stage is called the secondary oil recovery in which an external fluid is injected into the reservoir to support the depleted reservoir pressure from the primary recovery. The injected fluid can be either gas or water, and the widely used secondary oil recovery technique is waterflooding. After the secondary oil recovery process, a substantial volume of oil is still trapped in the pores of the reservoir rock and is economical to produce, hence leading to the tertiary oil recovery stage which is also known as enhanced oil recovery (EOR). A major focus in experimental efforts was directed toward understanding the role of geochemical reactions in CO 2 injection scenarios for such EOR schemes.
In EOR applications, when CO 2 is injected, the geochemical equilibrium is disturbed and the CO 2 can then flow in two phases: liquid and vapor. As the CO 2 contacts the formation water, it gets dissolved in the aqueous phase, while the rest of the CO 2 can either be dissolved in the hydrocarbon phase or roam freely as a vapor, Eq. 1. From a geochemical perspective, when CO 2 is dissolved in the aqueous phase, a dissolution process might occur, while ionic species such as HCO − 3 and CO 2− 3 are formed, which subsequently impacts the porosity evolution in the reservoir. In other words, CO 2 reacts with water to produce carbonic acid ( H 2 CO 3 ), which then dissociates into bicarbonate ( HCO − 3 ) and hydrogen ions ( H + ), Eq. 2. Generally, mineral interactions are dictated by the formation of the reservoir, where calcite is the main mineral in case of carbonate rocks. However, as the reservoir starts producing and the pressure decreases, the CO 2 in the gas phase will expand causing the CO 2 in the aqueous phase to exsolve. This causes the minerals to re-precipitate on the rock surface and might lead to significant changes in porosity and wettability of the porous media, and hence affecting hydrocarbon recovery. The set of equations below describe how gaseous CO 2 gas in the aqueous phase interacts with water to produce a number of different compounds: The schematic in Fig. 5 illustrates how pressurization affects the behavior of CO 2 in high pressure and low pressure environments. With pressure being low, the CO 2 might escape the gaseous phase and interact with the minerals in the aqueous phase causing their dissolution. As pressure is increased, CO 2 exsolves from the aqueous phase and the dissociated minerals begin to re-precipitate governed by the reversible reaction in Eq. 2. Other reactions might arise as well depending on the mineralogy of the reservoir rock and how it reacts with the present CO 2 . The injection of CO 2 in carbonate reservoirs causes the dissolution of the formation which greatly affects the capillary pressure and the relative permeability of the system. However, such changes have not been experimentally quantified and documented in well-established studies. The mechanism at which the intrinsic rock properties change is attributed to the enlargement of the pore throats. This creates new flow channels that enhances the connectivity of the system, as shown in Fig. 6.
The methods used to quantify and assess the interaction between the rock and carbon dioxide can be categorized into four groups:

Experimental studies
The experiments to investigate the rock-CO 2 interactions are performed on large and small scales and named total impact and specific experiments, respectively. In total impact experiments, the aim is to assess the immediate effects of CO 2 in the short term. For this purpose, CO 2 is typically injected to a host rock under controlled temperature and pressure, and the evolution of new species due the reactions in the system is observed. The work of [12] focused on evaluating the sealing efficiency 5 A schematic illustration of the pore and the reservoirscale levels of CO 2 mineral dissolution and precipitation with (a) and without pressure (b) [119] of a caprock upon CO 2 injections to handle the natural migration of CO 2 rich pore fluid through the caprock and in fault zones, while [22] studied a specific shale sample (Muderong Shale, Australia) with different clay minerals to obtain information on its sealing integrity. Many other similar experiments were also carried out and modeled in the laboratory through batch reactors [32,63,116,117] to understand the capture and geological CO 2 storage (CCGS) in deep aquifers and depleted oil and gas fields. On the other hand, specific experiments target selective minerals to be examined as opposed to the whole heterogeneous rock in the total impact studies. This approach is popular as it is scientifically more accurate when quantifying the effect of CO 2 on certain minerals. The results are highly advantageous in determining the reaction kinetics in the system and measuring the thermodynamic equilibrium constants associated with a specific type of mineral [33,86]. The experimental work of many scientists [9,70,91,113,114] provided great insight on how different minerals interact with the rock, but the developed results and models suffer from scalability issues to the field conditions. The measured rates of the kinetics of the reactions fail to represent field cases to date, as the research on this topic is still ongoing.

Natural analogues
The information obtained from natural CO 2 analogues could result in an better understating of CO 2 -water-rock interactions at the geological time scale [13,92]. Although such comparisons can provide important information about the interaction process, it is rather generic and qualitative. The studies can be broadly split into two types: studying the effect of CO 2 on the reservoir rock formation [16,50,118] or tracking the CO 2 source as it migrates throughout the subsurface [46,50,111]. However, it is almost impossible to replicate a geologic setting in different fields, as each location has its own set of unique properties and boundary conditions. This deems the analogue studies insufficient in providing critical analysis of the CO 2 role in rock-fluid interactions.

Field observations
The observation of already existing producing fields could provide insights on how CO 2 would react with the formation water and the reservoir rock. The data gathered about the rock properties, and the produced fluid is deeply investigated and analyzed to provide information on the potential behavior of CO 2 upon injection. Moreover, dedicated field-scale tests can be performed to track mineralogy changes and CO 2 migration in real time at different time scales. Pilot CO 2 injections have been carried out in many locations in the North Sea, Algeria, Canada [11,48,55,67,79], and the results were used to further support numerical modeling of reactive transport processes.

Numerical modeling
Since other methods come short when quantifying rock-fluid interactions at a large scale, numerical simulation is deemed essential. The models can be used to track the flow of CO 2 , investigate porosity evolution and quantify dissolution/precipitation of minerals in the porous Fig. 6 The mineral dissolution changes the pore structure by enlarging the throat diameters which is shown in the system to the right. As the pressure decreases, minerals start to deposit again on the rock surface [119] media. The main attraction in using numerical techniques for fluid-rock interaction scenarios is the capacity of the simulators to model sophisticated physical phenomena temporally and spatially [36], a feature that laboratory experiments lack in. Many studies have been performed to benchmark mathematical models to physical processes in the subsurface, specifically in the work of [30,61]. Multiple parameters were fine tuned in these studies such as the vertical permeability and the system boundary conditions, and the effects of CO 2 solubility in water were also considered [40,96].
The path to develop an efficient and robust geochemical model is not well paved. It requires a deep understating of the integrated processes and factors involved in the system under focus. Each scientist will look at the same system from his own background and point of view to determine the important features to be preserved by the conceptual model. But in order to understand how the coupling of reactions and transport is important to model a physical system, we first have to touch on some important concepts in reaction modeling.

What is chemical equilibrium?
When dealing with chemical reactions, we have to introduce three important concepts: chemical kinetics, equilibrium and chemical equilibrium. The chemical kinetics can be referred to as the science of studying and determining (both experimentally and theoretically) the rates of a reaction network. This reaction network can be formed of tens of simultaneous reactions where products and reactants are dependent on one another. The equilibrium condition is achieved when the forward and the reverse rates of each reaction equalize, and the reaction process stops. At this condition, chemical equilibrium is achieved where the net concentration of the components in the system is constant over time. In other words, we say that the system is in a minimum free Gibbs energy state at standard conditions. For a system at constant temperature and pressure, we can summarize the following results: 1. If G < 0 , the process occurs spontaneously.
2. If G = 0 , the system is at equilibrium. 3. If G > 0 , the process is not spontaneous as written but occurs spontaneously in the reverse direction.
The term G is defined as the change in the free Gibbs energy of the reaction and is denoted by: where R is the gas constant, T is the temperature, Q is the reaction quotient, and K is the equilibrium ratio. The terms Q and K will be discussed in detail in the section that follows.

How do we characterize a reaction?
Every reaction is essentially a reversible kinetic reaction. Thus, it is conceptually sound to model every reaction as a kinetic reaction. However, it might be extensive to gather experimental data about the rates of every possible reactions path in a system at different temperature conditions. Some systems might even prove challenging to establish its path and measure the rates in a laboratory setup due to the volatile nature of the components and their critical properties. Thus, we tend to put constraints in place that localize parts of the reaction system as we assume equilibrium conditions for some reactions. In reality, equilibrium is not achieved simultaneously, and the reversible reactions do not stop at equilibrium but rather the concentration change is constant. Let us assume a generic chemical reaction given by: At equilibrium, we define a quantity called the equilibrium ratio, in which equilibrium is achieved when the rates of the forward and the reverse reactions are equal. In the chemical reaction shown above, is written as: The ratio represents the concentration of the products to that of the reactants at equilibrium. This expression is conventionally called the law of mass action [109]. We illustrate schematically the significance of the equilibrium ratio of Eq. 6 for 3 different scenarios: The terms A 0 , B o and C o denote the initial concentration of each component in the 3 scenarios before any reaction occurs.
The illustration in Fig. 7 shows how the concentrations of the components change for the different cases of .
The ratio of change between the components is generally dependent on the stoichiometric ratios of the reactions, and the initial number of moles of the reactants. At initial conditions, the reactants are introduced and the reaction path (driven by the equilibrium ratio) favors either a forward or a backward reaction. For example, in case (a) of Fig. 7, the reaction is driven forward due to K eq >> 1 , and hence, A and B will fully react to form C with the backward reaction being of little significance. In case (c), where K eq = 1 , both products and reactants are present at equilibrium conditions and the ratios of the components depend on the stoichiometry of the reaction. When a reaction follows the law of mass action, we characterize it as an equilibrium reaction. These reactions are reversible and fast, as they occur spontaneously on a small time scale. This means that we can express them algebraically to relate with the reaction quotient as follows: where depends on the initial concentrations of the components in the mixture.
On the other hand, when the reaction has a strong dependency on time and is relatively slow and irreversible, we characterize it as a kinetic reaction. The term kinetic refers to the rate at which a chemical process occurs. If we rewrite the same generic reaction presented earlier in Eq. 6 as an irreversible reaction: we can then write the reaction rate as follows: The rate constant is denoted by k which is a reaction property that can be only determined experimentally and is highly dependent on temperature. Since a limited amount of rate measurements at different temperatures can done experimentally, we usually extend and interpolate the range for the rate constant using Arrhenius equation: Fig. 7 The concentration of the reactants and the products at initial and equilibrium conditions: a At high equilibrium ratios, the reactants are completely converted to C as the reaction moves forward b At low equilibrium ratios, the product is completely dissociated into A and B as the reaction moves reversibly c When the equilibrium ratio is close to unity, the reactions and the products tend to reach equilibrium at equal concentrations where E a is the activation energy (J/mol), R is the gas constant (J/k.mol), T is the absolute temperature, and A is the frequency factor (a constant to indicate the number of collisions with correct orientation that leads to products).
We can finally write the rate law with respect to A as: where x and y are the reactant order determined from the experiments. For simplicity, we can generally assume a reaction to be of first order (x = y = 0), but this treatment ultimately depends on the process being modeled. To summarize, reaction characterization is a complex process that requires rigorous experimental validation of the reaction mechanism. Since all chemical reactions can be, in essence, modeled as kinetic reactions, we selectively localize some of the reactions by assuming an equilibrium condition to further simplify the reaction network.

Kinetic or equilibrium?
It is up to the modeler to classify a set of reactions into kinetic and equilibrium reactions, as long as the integrity of the system is preserved. In theory, every reaction can be modeled as a kinetic reaction. However, it is almost impossible to do so because of the lack of experimental data to represent the kinetic rates of the reaction system. In this case, it is practical to assume some reactions are instantaneous and thus treat them as equilibrium reactions.
If we consider the reversible reaction of decomposition of N 2 O 4 (g) and dimerization of 2 NO(g), we write it as follows: In this reaction, equilibrium is not achieved instantaneously and concentrations stop changing at equilibrium but reactions do not stop. If we have a reaction vessel starting with 100% N 2 O 4 (g) that undergoes a full conversion to 2NO(g), this indicates that the kinetic rate of the dimerization reaction is zero. To test this out, we vary the kinetic reaction rate k for the decomposition reactions between 0.01 and 10 and we plot the concentration of the reactant and the product in Fig. 8. When the kinetic rate is high, the kinetic reactions approach the equilibrium point and thus can be treated as an equilibrium reaction.

Flow and transport in porous media
It is essential to understand the process of fluid flow in the subsurface geological structures, especially when multiple fluids are present as in the case of hydrocarbon reservoirs, see Fig. 9. This allows for proper modeling of the fluid path and preserving the integrity of the physical laws in critical processes such as underground water management, chemical weathering and oil extraction from deep reservoirs. Typically, transport is controlled by the phenomena Advection is the process where solutes move with the bulk fluid phase at the same velocity. Imagine you are next to a stream or a river, and you gently place a piece of paper on the surface of the water. The paper will start moving with the bulk phase (water) at the same velocity of the flow. Conversely, dispersion and diffusion are concentration driven and depict the extent of spreading of different components.
The generic continuity equation describing the concentration of a component with time for fluid flow in porous media, while accounting for the advective, diffusive and dispersive forces is represented by: where C i t is the accumulation term written as C i = j S j X ij , n p is the number of fluid phases, j is the density of phase j that depends on pressure, temperature and molar fraction, is the porosity , X ij is the mole fraction fluid species i in the fluid phase j, u j is the Darcy velocity for phase j, D ij is the dispersivity tensor for fluid species i in phase j, and q W i is the source/sink term of species i.
The Darcy velocity u j is given by Darcy's law: where k rj is the relative permeability of phase j, p j is the pressure of phase j, g is the acceleration of gravity, and K is the rock permeability. Typically, the Kozeny-Carman equation is used to dynamically represent the permeability in the continuum scale: where K 0 and 0 represent the initial permeability and porosity in the system, respectively, before any chemical reaction takes place. In case, the formation porosity is assumed to be constant throughout a certain physical process, and then, K = K 0

Coupling reactions with transport
The coupling of kinetic and equilibrium reactions with the mass balance equation of the fluid flow can result in a highly complex system of equations especially when compositional systems are modeled [4-6, 71, 123]. Since the reactive transport area is quite mature in the hydrology community, recent efforts in the work of [43-45, 62, 75, 95, 108] have been focused on utilizing the advances in the hydrogeological field to design a comprehensive framework for reactive flow in porous media for hydrocarbon reservoirs, see Fig. 10. In this section, we present the different variable formulations that are used to represent the continuity and reaction equations and highlight the key differences in each approach.

Natural variable formulation
The natural variable formulation was first developed by [31] allowing the implementation of compositional modeling in reservoir simulations with fully implicit schemes.
In this formulation, [31] consider n c hydrocarbon components plus the water phase which adds up to n c + 1 components, requiring similar number of mass balance equations to be solved. The primary variables are the pressure, phase compositions and saturations, while phase behavior relationships are used as local constraints. Many modifications and alterations to the natural variable formulation have been proposed to fit different physical applications.
In the following, we present different forms of the natural variable formulation for reactive transport with multiphase flow problems. The presentation and the analysis of the transport equations are based on the assumption that Darcy's flow is valid throughout the domain.

Governing equations of transport and chemistry
A component, in thermodynamics, is defined as a chemically independent constituent of a system. If we have a system that constitutes n s solid species and n f fluid species, then we can write the total number of components in the system as n c = n s + n f . Incorporating the reactive term to the mass balance relation in Eq. 14 results in: where N i is the accumulation term of component i, L i is the flux term accounting for the advective, diffusive and disp e r s i v e f o r c e s a n d i s w r i t t e n a s ( ∇ ⋅ ∑n p j=1 ( j X ij u j − j S j D ij ∇X ij ) + q W i ), v i,k is the stoichiometric coefficient of species i in kinetic reaction k, and v i,q is the stoichiometric coefficient of component i in equilibrium reaction q. The kinetic and equilibrium reaction rates are denoted by r k and r q , respectively, and n k and n q are the number of kinetic and equilibrium reactions, respectively. The total number of the reactions in the system can be then expressed by n r = n k + n q .
The accumulation term represents the change of mass of the fluid and solid species with time and is written as, where n jf is the number of fluid phases, j is the density of phase j that depends on pressure, temperature and molar fractions, is the porosity, X ij is the mole fraction of fluid (17) species i in the fluid phase j, s is the density of the solid phase, and X is is the mole fraction of the solid phase s. The natural variable formulation expresses the unknowns in the system in terms of the set To further simplify the system, we assume that the solid phase does not flow and phase partitioning only occurs when fluid components are present. This leads to a simplified set of equations based on Eq. 17: where C s is the concentration of the solid components s, n f is the number of fluid components, and n s is the number of solid components ( recall n c = n f + n s ).
Subsequently, the natural variable set of unknowns is then rewritten as, The governing equation represented in Eq. 17 is the basis of all the simulators discussed earlier in the introduction. The main differences only lie in the treatment of the kinetic and equilibrium reactions. If all the reactions in the system are kinetic reactions, the rate of the reaction is then expressed in terms of the natural variable set and can be solved with the flux term or replaced with a rate law equation similar to the ones discussed in section 4.2. On the other hand, equilibrium reactions are more complicated due to their fast nature and are usually expressed by algebraic equations. This means that each equilibrium reaction will add additional unknowns to the system, and computational complexity overall [74]. In the hydrologic simulators, the number of equilibrium reactions is reduced by formulating the problem in terms of components. A component in this case is defined as the combination of species that are not affected by the equilibrium reactions in the system and is different from the definition of a component in reservoir simulations [80]. This makes such a concept difficult to apply in multiphase flow context, and hence, an alternative approach is needed to both simplify the system and preserve the characteristics of the reactions.

Element-based formulation
Coupling geochemical reactions with transport equations for hydrocarbon production applications poses v s,q r q (s = 1, ..., n s ) complexities on the solved system due to the differences in the time scales of the transport and chemistry [110]. Such differences would typically result in a stiff system of equations, and hence, an effective strategy needs to be adopted to solve or simplify such systems. The work of [23,106] discussed how different time scales could affect the quantification of the deposition process in the porous media at the pore, core and field level using analytical and numerical modeling techniques. The results showed the influence of the Damkohler number (Da), a dimensionless quantity that relates the chemical reaction time scale to the transport phenomena rate occurring in a system, on the properties of the system. If D a > 1 , reaction is dominant and clogging effects are exacerbated and evident in the porous media as the variations in the concentration field are significant. However, and for D a < 1 , convective forces are dominant and the porous medium properties change smoothly. At a field level, the prediction of simulations can be improved by using relevant core-scale transport properties of reactive flow for each rock-type and reactive flow regime [106].
As mentioned before, one could linearly transform the transport equations to eliminate the dependency on equilibrium reactions (and hence the fast nature of the reactions) using an element-based mass balance equations [20,99]. The governing equations 20 and 21 can be further simplified and written in a concise matrix format: where is the vector of total component concentration, = (N 1 , ..., N nf , C1, ..., C n s ) T ; is the flux vector, = (L 1 , ..., L nf , 0, ..., 0) T ; is the reaction term vector, = r k r q and is the stoichiometric matrix that holds the coefficients of the components in each reaction with dimensions n t × n r , = S k S q where the stoichiometric submatrix S k corresponds to kinetic reactions and S q to equilibrium reactions. The total number of reactions in the system, n r , consists of kinetic n k and equilibrium n q reactions of which we can safely assume that they are linearly independent and balanced ( n r = n k + n q ) [43]. By definition, an element is an atom or a compound that does not partition into smaller entities by chemical reactions [44]. Let us take a matrix n e ×n c that holds the stoichiometric values for each element in each component. When the set of reactions is assumed to be linearly dependent, then the number of components is equal to the number of reactions and the number of elements in the system ( n c = n e + n r ). Thus, if we multiply by the component stoichiometric coefficient matrix , since both matrices are nonzero matrices and do not have a full rank, we get: Let us illustrate the relationship above with an example of a set of independent balanced reactions below (the reaction can be either kinetic or equilibrium) : The elements Thus, we managed to eliminate the reaction term from the mass balance equations (whether they are kinetic or equilibrium) which adds an advantage to the element-based formulation. We end up with a system of conservation equations for each element involved in the reactions. In other words, we solve for Eq. 27 plus n r reaction relations instead of solving Eq. 23 directly.

Reduction of reaction rates
The approach, in theory, is similar to the element-based formulation reduction procedure, as we have to linearly transform the system to eliminate the equilibrium rates but through multiplying the element stoichiometric matrix, with the canonical or reduced row echelon form of the stoichiometric coefficient matrix, denoted now as ̄ . The applied mathematical manipulations to obtain the canonical matrix depend on the kinetic and equilibrium reactions in the system, and the choice of primary and secondary components. Different forms of ̄ for specific reaction sets have been discussed in the work of [45,80]. The product of ̄ would then result in where is an identity matrix. If we rewrite Eq. 23 in terms of ̄ and multiply by the matrix , we get: The resultant system signifies that each kinetic reaction rate only appears in one equation due to the linear (24) n e ×n c × n c ×n r = 0 n e ×n r (25) combination of the component conservation equations.
But since the kinetic and equilibrium reaction relations are provided separately from the mass balance equations, one should have the knowledge to model different types of chemical reactions and their associated rates. Next, we present a generic way to treat chemical reactions despite the fact that such treatment is usually simplified depending on the problem at hand.

Generic reaction implementation
In order to implement the reactions into a Jacobian, we have to obtain the derivatives of the reaction term with the set of variables . It is essential to denote that the reaction term r k depends on the concentration/activity of the fluid species i in phase j and the concentration/ activity of the mineral solid phase s ( r k (C ij , C s )). The reaction is then characterized in terms of both fluid and solid concentrations/activities.

Concentration-based reactions
If we take the derivatives of the reaction term with the respect to the set of unknowns , we obtain the generic derivatives for kinetic and equilibrium relations, respectively: In order to derive the generic reaction relations with respect to the concentration, one must assume a form of the reaction rate, see Table 1, for example.
Let us assume that we have kinetic reaction that follows the law mass of action. In this case, Eq. 29 must be used to obtain the derivatives and the reaction rate would be in the form of: where K k is the reaction rate constant, m and n are the stoichiometric coefficients of the element, and the derivatives of r k with respect to for Eq. 29 are : We also relate C ij to set by: The derivatives of the fluid concentration with respect to the unknown variable set will then be: where and are target phase and species derivative indices and is the Kronecker delta function.
Similarly, we calculate the derivatives of C s with respect to the variable set. Since we have a solid phase, only the derivative of C s with itself is nonzero, where is a mineral component index.
If Eqs. 32-33 and 35-37 are plugged back to Eq. 29, the resultant expression can then form the basis of the Jacobian representing the reactions in the system. If the equilibrium reactions are present in the system, a similar derivation procedure is adapted through Eq. 30 to obtain the required derivatives for each equilibrium reaction rate, Q q .

Activity-based reactions
In some problems, the solubility of components in the gas phase is better represented using an activity approach, specifically in precipitation/dissolution processes. Thus, the generic reaction equations must be expressed in terms of activity coefficients. We denote the activity coefficients of fluid species i in phase j and the solid species s by a ij and a s , respectively. The generic derivatives for kinetic and equilibrium relations, respectively, are written as: r k a s a s x Since the activities of the solid phase are set to unity, then the derivative of a s with respect to is 0. The kinetic rate is represented in terms of the activity of the components in the fluid phase. An example rate law for such processes is written as: where A k is the reactive surface area used by reaction k, K k is the rate constant, and K eq,k is the equilibrium constant. The activity relations can be expressed in two distinct forms based on the state of the fluid, as fugacity is used when the components are dissolved in the gaseous phase, while molality is used when the aqueous phase is the solvent. The fugacity is defined as an empirical quantity to account for the fluid deviation from ideality and is measured by the effective pressure of a gas to express its tendency to escape or expand. The dimensionless activity relations then become, where f ij is the fugacity of species i in phase j = g , P 0 is standard pressure, X ij is the mole fraction of solute species i in phase j = w , X jj is the molar fraction of water in the aqueous phase ( j = w ), and ij is the activity coefficient for species i in the aqueous phase j = w.n. The constant value 55.508 represents the number of mole of water per kilogram, and the activity coefficient ij is calculated from empirically correlated models similar to [57,85].
Based on the presented definition of activity in Eqs. 42 and 43 and the rate equation in 41 , we can then obtain the derivatives of these expressions with respect to the variable set in a similar fashion to Eqs. 32-33 and 35-37. The derivatives are substituted in Eqs. 39 and 40 to obtain the final form of the Jacobian equations.
a ij =55.508 ij X ij X jj , j = w, i ≠ w (aqueous phase)

Summary of equations
In summary, we have 2n p + n s + n p n f variables requiring n e + n r + (n p − 1)n f + 2n p equations. We recall that n r is the number of chemical reaction relations n r = n k + n q , n p is the number of phases, n f is the number of fluid components and n s is the number of solid components where n c = n s + n f . When a system is fully characterized, this implies that n c = n f + n s = n e + n r , and thus, the total number of equations reduces to be equal to the total number of unknowns. We would like to also point out that the system is closed with a set of constraints as follows: In order to simplify the final system, it is possible to split the governing equations, reaction relations and local constraints into primary and secondary equations as we show in the following section.

System reduction: selection of primary and secondary equations
We first solve for the primary unknowns, which are then substituted into the secondary set to compute the secondary variables. This procedure does not alter the total number of the governing equations, but rather selectively categorize the equations into different sets to increase the performance of the simulator at the liner solver stage. In order to decompose our governing equations into two categories, [44] suggested two approaches: reactions coupled and local reactions decoupled. In the reaction coupled approach, we consider our mass balance equation including the reaction term (both kinetic and equilibrium) to be our primary system, while the local constraints are treated as secondary. The primary variables in this selection would be the pressure (P), phase saturation ( S j ) and some phase composition ( X ij ), while the rest of ( X ij ) and the solid-phase composition ( X is ) are the secondary variables. Similarly, for the second approach when local reactions are decoupled, the primary and secondary variables selection is the same except for the number of the primary equations. The equilibrium reactions are separated from the primary set and treated as secondary equations.
In some cases, further reductions can be performed by assuming localized chemical reactions, but this solely depends on the reaction system under study. Table 2 summarizes the equation selection criteria discussed earlier.

Solution procedure
The Jacobian is constructed based on the selection of the primary and secondary equations and is solved in a fully implicit scheme. Although the fully implicit approach guarantees a stable solution, the strong coupling of different physical phenomena in the system results in a stiff matrix. Moreover, this approach might be memory intensive and pose a computational complexity in the simulator. This elevates the need to develop unique pre-conditioners for the system based on the modeled physical processes and the numerical discretization scheme that is used [81]. In a reactive transport system, the primary equations are solved at each Newton step which necessitates the computation of the derivatives of the governing reaction with respect to the variable set at each iteration. Thus, decoupling might prove handy in simplifying the procedure of obtaining such derivatives and speeding up the computation time. After the primary variables are calculated, we move to perform thermodynamic analysis of the system to study the phase stability and determine whether phase disappearance/reappearance will occur at the next Newton iteration. Finally, the secondary variables are then computed. The procedure presented in algorithm 1 summarizes the solution steps.
After we have fully presented the natural variable formulation and the different approaches in representing the governing equations and the reactions relations of a system, we introduce a molar variables basis for the system formulation. In the following, we will discuss two different molar variables formulations and highlight the key differences of each approach.

Molar variable formulation
The molar variable formulation was first introduced by [97] and fully tested for compositional reservoir simulations by [62]. First, we define a new porosity term, ̂ , that refers to the reactive porosity defined as the volume occupied by the solids that do not participate in chemical reactions and related to by: Reactions coupled Mass conservation equations [28] Phase equilibrium and local constraints [44,45,46] Reaction relations [29,30] Local reactions decoupled Mass conservation equations [28] Equilibrium reaction relations (localized) [30] Phase equilibrium and local constraints [44,45,46] Kinetic reaction relations [29] Kinetic reaction relations [29]

Overall-composition formulation
Now, we would like to express the governing equations in terms of overall-composition variables using the following set: where z c is the overall-composition fraction, X ij is phase composition, and j is the phase molar fraction defined as: The fluid saturation S j is expressed as: We relate the overall composition of each component to the molar fraction of the fluid by: We also define T as the total component density: The above relationships allow us to propose and alternative form of N i in Eq. 18 as follows, This form allows us to write the governing equations similarly to Eq. 28 and then reduce the n c component conservation equations into n e element conservation equations using the matrix , .., n c &j = 1, ..., n p ) The reactions are treated independently similar to the procedure outlined in Sect. 6.1.4. Table 3 summarizes the primary and secondary variables of the system. The total number of variables in the system is n c n p + n c + n p + 1 closed by the element conservation equations, reaction equations and local constrains which will be discussed next. Note that the number of variables is equal to the number of governing equations. The local constrains are defined as follows: (1) Phase Equilibrium Note that not all z i values are primary but rather those related to n e + n k − 1 equations. If we further decouple the chemical reactions to include the equilibrium reactions in the secondary equations, we end up with a set that includes the residuals of the local constraints, phase behavior and equilibrium reactions equations. The solution procedure is quite different than the one discussed before for natural variables. Here, we assume that the primary variables are known (P and part of z i ), and the goal is to determine j , X ij and the rest of z i from the secondary set of equations. The procedure involves two steps; the first step is to determine the number of existing phases in a given cell. This follows with flash calculations to determine the composition of the existing phases ( X ij ) as the second step.
After the secondary variables are obtained from the flash calculations, we update the residual of the primary equations to check for convergence. If the results are unsatisfactory, then we are required to calculate the partial derivative of the secondary variables with respect to the primary variables: This derivative is computed using the results from the last converged step , where then the Jacobian of the primary equations is reconstructed and solved for the primary variables. This process is repeated until Newton iterations converge. Algorithm 2 shows the steps of the process.
The main advantage of this formulation is that the number of equations and variables is the same in any given cell. However, phase computations are required for each cell at every step, and thus increasing the number of flash calculation iterations.

Element molar formulation
In this formulation, we take the elements in the system as the basis of the governing equations as opposed to components. The number of variables is similar to the overallcomposition formulation with additional n e variables. The following variable set is used: where z e is the ratio of moles of an element to the total moles of an element.
The element molar formulation can be easily expressed using a set of relations to perform component-to-element transformation from z i to z e : where E T is the total element density, equivalent to the total moles of elements divided by pore volume. The total element density is expressed as: where E ei is the number of element e in component i. We apply these changes to the accumulation term of Eq. 54 to write the governing equations in terms of the element molar formulation by dividing with E T : The solution procedure is quite similar to the overall-composition formulation. Here, we assume that the primary variables are known (P, z e and part of z i ), denoted as the primary set e,p . Table 4 summarizes the primary and secondary variables of the system. The goal is to determine j , X ij and the rest of z i from the secondary set of equations. In this formulation, the last z e value is not solved for using the primary equations but rather computed from the constraint: When we formulate the secondary equation set , we add an additional constraint to represent the component-toelement transformation in the system. Since phase behavior computations are performed on components, we need to constantly combine the elements to form components when needed. P, z e (e = 1, … , n e − 1) , z i (i = 1, … , n k ) z i (i = n k + 1, … , n c ), X ij , j In molar formulations, phase behavior computations are straightforward if chemical reactions are not involved. Typically, the pressure and the molar composition z i are known throughout the calculations, and hence, a twostage procedure is easily implemented. First, we perform a stability test to determine the number of phases in a grid block [49,78], followed by traditional flash or negative flash calculations to determine the molar fraction of each component in each phase. However, when chemical reactions are implemented, the molar composition of the mixtures is not totally known. In particular, n e + n k − 1 components of z i are known, which might limit phase behavior procedures in existing reservoir simulators. The chemical reactions pose an additional constraint on the parameters represented by chemical equilibrium, a term that describes the interactions of multicomponents in multiple phase. Adding this constraint to the existing phase equilibrium requirement presents critical challenges to phase behavior computations. In brief, chemical reactions need a special combined thermodynamic and chemical equilibrium approach to ensure phase computations are accurate and representative of the physical system being investigated.

Which formulation is better: molar or natural?
Traditionally, natural variables have been widely used in commercial reservoir simulators due to the ease of implementation and the simple phase behavior algorithms. However, as the more complex physics are required to be modeled, molar formulation started to be implemented for specific simulation scenarios. Each formulation has its own advantages and drawbacks, which we summarize in four main points.

Jacobian construction
Typically, the derivatives of the governing equations are required to construct the Jacobian matrix and solve for the primary unknowns. When using natural formulation, computing the derivatives is a straightforward process since the governing equations are naturally expressed in the terms of the set of the unknowns (hence the name of the formulation). However, in the molar formulation, the expression of the derivatives is much more complicated as we have seen in Eq.59. In a reservoir simulator, such approach would increase the processing time of setting up the Jacobian matrix and thus reduces its desirability.

Variable substitution
The size of the Jacobian matrix is decided by the number of equations and unknowns in the system. Thus, the Jacobian has to dynamically represent the system when any changes occur. When natural variables are used, the primary variables set has to represent the phases present in the system. As the simulation progresses, phase appearance/disappearance is experienced in some cells, and thus, the number of unknowns varies accordingly. Such behavior adds unwanted complexities in a reservoir simulation when performing phase behavior calculations, especially when parallel processors are utilized. The partitioning of the Jacobian into different computer cores is not straightforward as the size of the system varies with each time step. On the other hand, molar formulation does not suffer from phase appearance/disappearance issues since overall variables are used. This means that the number of unknowns is the same in each grid at any given time step. The phase behavior is decoupled from the primary set and is solved separately.

Flash calculations
In the natural variable formulation, the phase behavior is solved together with the governing equations to determine the phase equilibrium conditions. However, in the molar formulation, the flash calculations are performed separately and referred to as exact flash. This treatment ensures that accurate depiction of the phase state of the system is achieved as the time progresses in a reservoir simulation. However, this approach has to be further investigated when chemical reactions are involved.

Computational cost
The phase behavior computations in the molar formulation increase the complexity of the system and subsequently increase the cost. However, this should not be the sole factor in picking a method over the other. Factors such as the time-stepping size and Newton iterations vary with each formulation based on the system being studied. Furthermore, decoupling chemical equilibrium reactions in the natural formulations will reduce the computation cost further and this should be taken into consideration when assessing the performance of both formulations. The advantages and the disadvantages of the natural and molar variable formulation are summarized in Table 5.
Earlier, we presented the advantages and disadvantages of utilizing molar and natural formulations for coupled reaction systems. Ultimately, the choice of the method will be based on the physics of the problem and the desired efficiency of the simulation code. Next, we present two computational cases from the literature that relies on the previously discussed formulations, and we briefly comment on the results of each numerical case.

In situ upgrading of shale oil
In the first case, we show an in situ upgrading of oil shale model by [43]. The process is highly dependent on temperature which directly affects the rates of the kinetic reactions as compositional flow is experienced in the reservoir. [43] modeled the problem using a component-based formulation that we presented in Sect. 6.1.1 using natural variables. In the model, the kinetic chemical reactions are strongly coupled with the mass and energy balance since they are highly dependent on temperature, while average values are used to represent the porosity and permeability changes. Furthermore, the decomposition of kerogen into multiple mixtures yields four different phases (gas, water, oil and solid prechar) in which the solid phase is assumed to not occupy any pore space. The decomposition reaction of the kerogen adapted from [64] and reported in [112] is as follows, where the symbol IC indicates lumped components. All the reactions are assumed to follow the elementary rate law in Eq.31, and the rate constant K k is given by Arrhenius law in Eq. 11.
Moreover, since most in situ upgrading processes occur above the boiling point of water, the problem is reduced to a gas-oil equilibrium analysis. However, if the water production is significant, three phase flash calculations are required to thermodynamically represent the system. The unknowns to be computed in the system are pressure, temperature, mole fractions of each fluid component and solid concentrations, while localized reactions are strongly coupled. The model of [43] was used to validate the results of the Mahogany Demonstration Project-South (MDP-S) pilot in situ upgrading project [47] where the pore space of the matrix is initially assumed to be saturated with brine or gases. Heaters are spaced around the producers to facilitate the decomposition of kerogen where oil starts accumulating in the pore space and gets produced by the well, see Fig. 11.
The results of the simulation are presented in Fig. 12 where high initial gas production is noticed due to the heating effect to accelerate the decomposition of the kerogen into smaller carbon chains. The plots of the cumulative production show a great agreement with the field results. The importance of using a reactive transport approach in this problem is related to the decomposition of the kerogen into smaller hydrocarbons that need to be tracked, while the bulk fluid phase flows in the reservoir. Since the water phase is not of great importance in the in situ upgrading, the component-based formulation (presented in Sect. 6.1.1) is used to monitor different alkanes in the system as the temperature changes and the kerogen is consumed.

Coupling high-order linearization schemes with molar variables for reactive transport modeling
In the previous sections, we have extensively detailed the different formulations used to set up a reactive transport problem mathematically and discussed the pros and cons of each approach. The selection of the formulation that optimizes the problem is essential as the main objective of this procedure is to drive the residuals of the secondary equations to zero at each Newton iteration accurately and efficiently. However, this process is quite expensive especially when complicated derivatives are computed as in the molar formulation. Fig. 11 The heater spacing around the producers is shown in a while the gridding of the domain is depicted in b (reported in [43])

Fig. 12
The plots show oil and gas production rates (left) and the cumulative oil and total production (right) for the simulation of in situ upgrading of oil shale in the MDPS-S project [43].
The numerical and field production data are in good agreement, providing the possibility of modeling and designing in situ shale upgrading projects using the component-based approach To overcome this issue, [66] developed a multilinear state-dependent operators based on an operator-based linearization (OBL) technique to simplify the computation of the Jacobian derivatives. The technique is based on dividing the grid into smaller subdomains, encompassing one or more grid cells, and calculating the physical properties at specific points in the subdomains. This results in pre-computed properties stored in specific operators which are used to compute the derivatives inside the Jacobian. The derivatives are easily computed at any location in the grid through interpolating the values at the populated points, see Fig. 13 In OBL, the governing equation (Eq. 63) is discretized and rewritten in terms of state-based (physical properties of fluid and rock) and space-based (properties transformed in space) operators: where is a state-dependent parameter, is a spacedependent parameter, and and n are nonlinear unknowns at the current and the previous time step, respectively.. The operators for a simple equilibrium reaction case are given by, where c r is defined as the rock compressibility, T ab is the geometric part of transmissibility (which involves permeability and the geometry of the control volume) , and e ic is the amount of element i in component c.
The values of the operators and can be computed beforehand and prepared in tables at different compositions for a specified pressure range. It is important to note that the summation of z i of all elements equals to unity, and thus, one can produce tables of the operator values at a given resolution as a pre-processing stage. It is also possible to integrate an adaptive OBL scheme within the process of the simulation adding to the efficiency of the reservoir simulator. This technique is quite useful in reducing the complications of the phase behavior calculations and has been rigorously tested for complex geothermal and multicomponent flow and transport problems with buoyancy [65], and for reactive flow of species in reservoirs [62]. The working flow of the OBL approach is shown in Fig. 14.
The work of [62] presented a 1-D case of calcite dissolution using element formulation and molar variables. The dissolution and precipitation reactions were given by, The test was implemented using the OBL approach, and the equilibrium constants for the equilibrium reactions were scaled up to show the effects of the reactions on the overall system. The initial elements present in the reservoir are water and some other ions, while CO 2 is injected continuously. The simulation results highlight the effect of CO 2 on the behavior of ions in the reservoir, and the variations in their concentrations after 2000 days of simulation, as shown in Fig. 15.
The injection of CO 2 starts at distance 0 m causing the vaporization of water which in turn accelerates the calcite deposition. This causes an increase in the reactive porosity The domain of the reservoir grid. The black points represent the locations at which the OBL operators are computed. The red points A and B are interpolated based on the pre-computed operators at the surrounding black points, thus easily determining the value of the derivatives at the red points of the system and in turn impacts the permeability of the reservoir. The production (precipitation) and dissolution of ions will contribute to the transport of fluids in the system as more and more CO 2 is injected. The molar formulation was used in the latest case opposed to the approach of [44] to easily couple it with the OBL scheme and significantly improve the run-time of the simulation. This also shows that the molar framework is capable of handling multiple equilibrium reactions when integrated with high linearization schemes such as OBL.

Closing arguments and the way forward
The continuous increase for the need of hydrocarbons in the marketplace has pushed for constant improvements on the performance of the existing petroleum reservoir simulators physically and numerically. The area of reactive transport where chemical reactions greatly affect the flow of different fluids in the reservoir and the characterization of rock-fluid interactions through precipitation and dissolution processes are a new attraction for current research in the petroleum applications of reactive transport, which In this work, we described the physical significance of geochemical equilibrium and the reaction systems on the surrounding environment in porous media. A detailed procedure was outlined on how to mathematically implement the reactions in the framework of an existing reservoir simulator for both natural and molar variables formulations. A brief comparison was drawn on the performance of the two formulations and was supported with a review of two numerical cases from the literature. Ultimately, the choice of the formulation is dictated by the physical problem to be modeled where a trade-off has to be made between accuracy and efficiency.
In the future, we aim to integrate these reaction implementation techniques with a mimetic finite difference (MFD) discretization scheme [5] and a parameterized set of governing equations using the OBL method while using unstructured grids and full tensor permeability [2]. Upon the completion of such a simulator, it would be able to simulate field problems of hydrocarbon recovery with realistic scenarios such as heterogeneous 3-D aquifer and artificial groundwater recharge, CO 2 storage and sequestration and in situ upgrading of shale oil.
Declarations The authors have no conflicts of interest to declare that are relevant to the content of this article Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.