IORSim: A Mathematical Workflow for Field‑Scale Geochemistry Simulations in Porous Media

Reservoir modeling consists of two key components: the reproduction of the historical performance and the prediction of the future reservoir performance. Industry-standard reservoir simulators must run fast on enormous and possibly unstructured grids while yet guaranteeing a reasonable representation of physical and chemical processes. However, computational demands limit simulators in capturing involved physical and geochemical mechanisms, especially when chemical reactions interfere with reservoir flow. This paper presents a mathematical workflow, implemented in IORSim , that makes it possible to add geochemical calculations to porous media flow simulators without access to the source code of the original host simulator. An industry-standard reservoir simulator calculates velocity fields of the fluid phases (e.g., water, oil, and gas), while IORSim calculates the transport and reaction of geochemical components. Depending on the simulation mode, the geochemical solver estimates updated relative and/or capillary pressure curves to modify the global fluid flow. As one of the key innovations of the coupling mechanism, IORSim uses a sorting algorithm to permute the grid cells along flow directions. Instead of solving an over-dimensionalized global matrix calling a Newton–Raphson solver, the geochemical software tool treats the species balance as a set of local nonlinear problems. Moreover, IORSim applies basis swapping and splay tree techniques to accelerate geochemical computations in complex full-field reservoir models. The presented work introduces the mathematical IORSim concept, verifies the chemical species advection, and demonstrates the IORSim computation efficiency. After validating the geochemical solver against reference


Introduction
Reservoir simulators are designed and optimized for modeling fluid flow in complex geological structures with decades of injection and production histories.The industry-standard software programs allow a reliable computation of physical parameters such as phase pressures, saturations, and temperature to test field development scenarios, drainage strategies, and recovery optimization.
Rock-fluid interactions influencing the subsurface flow have mostly been overlooked, not because they are irrelevant but because geochemical reactions introduce numerical complexity that significantly affects the computation efficiency.Close to wellbores, rock-fluid interactions lead to scaling, and specialized software is used to monitor and plan scale treatments (Sorbie et al. 1991).However, rock-fluid interactions are not limited to the wellbore vicinity and may occur deep inside geological formations.
Field operators explore and execute improved and enhanced oil recovery methods as high recovery factors characterize most hydrocarbon reservoirs.Tertiary recovery methods include the injection of miscible and immiscible fluids, thermal applications, or substances that cause chemical reactions to increase oil recovery (Lake et al. 1989;Muggeridge et al. 2014).Moreover, underground gas storage is attracting increasing attention in the context of energy transition.The injection of carbon dioxide and hydrogen into porous underground media leads to fluid-fluid, fluid-rock, and/or microbiological interactions that significantly impact storage operations (Kelemen et al. 2019;Strobel et al. 2020).
Seawater is commonly used as a cost-effective fluid to displace oil in the Norwegian continental shelf (NCS).Although primarily injected for pressure maintenance, seawater may interact with the reservoir rock, induce formation compaction, or, ideally, improve oil recovery over waterflooding.To identify the impact of brine composition on oil recovery, experimental studies initially set a benchmark by evaluating the recovery efficiency of conventional waterfloodings (Strand et al. 2006).A modified brine then replaces the imbibing/injection fluid to estimate improved/enhanced oil recovery potential.Although many successful experimental enhanced oil recovery studies have been reported, upscaling the experimental findings to field-scale is complex.As a consequence, the number of implemented improved/enhanced oil recovery field operations remains limited (Bartels et al. 2019).
The Ekofisk field is one of the largest oil fields on the NCS and has exhibited significant compaction during more than 50 years of production.Apart from pressure depletion, Doornhof et al. (2006) attributed one-quarter of the Ekofisk compaction to rock-fluid interactions.Experimental studies have shown that selecting injection water compositions can affect compaction rates (Madland et al. 2011).Moreover, as the Ekofisk field consists of several distinct formation brines, monitoring the water ion concentrations provides information on aquifer drive, interconnection of formations, and well communication patterns.

Geochemical Modeling
In addition to solving the flow equation, reactive transport simulators must address nonlinear chemical species balances.Generally, two main formulations govern chemical equations in reactive transport modeling (Anderson and Crerar 1993;Bethke 2022): 1. Direct minimization of the Gibbs free energy function.2. Formulating a set of nonlinear equilibrium reactions with the law of mass action.
Both formulations lead to a highly coupled, nonlinear set of discrete equations that sequential or global implicit approaches can numerically solve (Yeh and Tripathi 1989;Saaltink et al. 2001;Steefel and MacQuarrie 2018).Sequential approaches use operator splitting to decouple the transport equations from the chemical reactions (Kirkner and Reeves 1988; Reeves and Kirkner 1988).On the one hand, decoupling the transport and chemical reaction calculations results in splitting errors, which constrains simulations by putting an upper limit on the time steps (Valocchi and Malmstead 1992).On the other hand, sequential approaches are simpler to implement.Simulators that use a sequential approach are Hydrogeochem (Yeh et al. 2012), THOUGHREACT (Xu et al. 2008), and the onedimensional transport solver in PHREEQC (Charlton and Parkhurst 2011;Parkhurst and Appelo 2013).
Unlike sequential reactive transport reservoir simulators, global implicit numerical schemes solve the transport equation and chemical reactions simultaneously.Characterized by good stability and robustness, global implicit approaches can address large time steps without encountering convergence problems.Although global implicit schemes avoid operator splitting errors, handling reactive multiphase and multicomponent systems results in tremendous memory demands.Yeh and Tripathi (1989) concluded that a reasonable utilization of global implicit approaches is restricted to one-dimensional simulation problems.Although more advanced global implicit approaches have been developed recently, the number of grid blocks processed remains significantly below the number of grid cells typically used in reservoir engineering applications (Saaltink et al. 2001;Kräutle and Knabner 2007).Examples of simulators that use the global implicit method are MIN3P (Su et al. 2021) and PFLOTRAN (Hammond et al. 2014).
Allowing the computation of large time steps, reservoir simulators employ implicit schemes to calculate unknown phase pressures and saturations.In the case of compositional models, implicit methods apply flash calculations to describe the partitioning of the aqueous, oleic, and gaseous phases into individual components (Trangenstein and Bell 1989).Typically, the Newton-Raphson method is used to iterate toward a solution by solving all unknowns simultaneously.While implicit approaches are feasible for black-oil models, the processing of multiphase and multi-component systems tends to introduce numerical instabilities and high computation costs (Rasmussen et al. 2021).
The literature lists IPhreeqc-UTCHEM, Reaktoro-Firedrake, and TOUGHREACT, among others, as examples of geochemical solvers that have been successfully coupled with flow simulators.IPhreeqc provides an interface that allows external programs to communicate with the PHREEQC geochemical packages.A detailed description of the IPhreeqc and UTCHEM coupling capabilities is presented in Korrani et al. (2015).To reduce computation time, IPhreeqc performs chemical equilibrium calculations only at the initial time step and subsequently adjusts the solution's speciation without full recalculations.To maintain material balance, IPhreeqc has to track the charge imbalances and the transport of hydrogen (H), oxygen (O), and other elements throughout the reactive transport computations (Korrani et al. 2015).Reaktoro is an open-source computational framework that employs numerical methods to calculate various reactive processes governed by chemical equilibrium and/or chemical kinetics (Leal 2015).In addition, Reaktoro offers on-demand machine learning methods to reduce computation times (Leal et al. 2020).The literature provides examples of Reaktoro coupled with Firedrake and DuMuX to simulate geochemical processes in porous media (Kyas et al. 2022;Koch et al. 2021).TOUGHRE-ACT is a multiphase, non-isothermal reactive transport simulator using shared memory parallelization.The software enhances the capabilities of TOUGH2, the initially developed multiphase flow solver, by incorporating a range of geochemical reactions such as aqueous reactions, mineral dissolution/precipitation, and gas reactions (Xu et al. 2006).

Adding Geochemistry to Existing Reservoir Models
Most hydrocarbon reservoirs on the NCS are modeled and history-matched using SLB Eclipse (Schlumberger 2002).While providing simplified and pre-defined formulations to simulate selected geochemical problems, the software is not designed to handle field-scale reactive transport problems.IPhreeqc-UTCHEM, Reaktoro-Firedrake, and TOUGHRE-ACT, among others, provide geochemical formulations; however, lack the ability to integrate realistic and history-matched reservoir models.This paper introduces IORSim, an add-on software to existing reservoir simulators, to address two main challenges of reactive transport simulations: 1. Providing a numerically robust and computationally efficient approach for solving fieldscale reactive transport problems.2. Allowing a seamless computation and visualization of geochemical reactions in existing reservoir models.
After computing the velocity fields of the fluid phases (e.g., water, oil, and gas) in a standard reservoir simulator, IORSim employs topological sorting to identify the upstream and downstream neighbors of a grid block (cf.Sect.2.1).Instead of solving a three-dimensional nonlinear IORSim species balance simultaneously, the identification of upstream and downstream cells allows IORSim to decouple the geochemical reaction calculations into an iterative sequence of local sub-problems.Solving the nonlinear geochemical reactions as a set of one-dimensional problems improves numerical stability and reduces computation times.
Besides the decoupled geochemical reaction species computations, IORSim offers additional numerical tools to tune and improve the reactive transport modeling.Depending on the simulation objectives, the software provides two specialized coupling mechanisms (Sect.2).First, the unidirectional high-speed forward simulation mode allows the post-advection of chemical components to visualize species distribution, water chemistry, and mineral interactions.Secondly, if the geochemical reactions affect the fluid flow, the bidirectional IORSim backward coupling calculates relative permeability curves to adjust the global fluid flow at each time step.If a large number of chemical components and geochemical reactions interfere with each other, IORSim provides basis swapping and splay tree techniques to reduce computation times (Sect.3.5).In a typical fieldscale reservoir model, large grid blocks tend to introduce numerical dispersion and distort geochemical effects.Instead of reconstructing the entire geological model, IORSim allows the definition of one or several finer grid resolution regions to improve the geochemical species balance calculations (Sect.3.3).
Reservoir models are typically optimized over years and include long injection/production histories, fluid properties, and well data within complex geological structures.Moreover, assisted history-matching software tools are commonly used to accurately represent the subsurface flow.Thus, depicting existing reservoir models in external reactive flow simulators requires substantial effort.By integrating IORSim into existing reservoir simulators, IORSim allows a seamless and straightforward representation of geochemical subsurface phenomena.
At present, IORSim features a water diversion module (sodium silicate injection), partitioning, non-partitioning tracer module, and a sophisticated geochemical solver that accounts for mineralization, ion exchange, surface complexation, and other chemical reactions (cf.Sect.4.1).
In a previous study, IORSim was used to correlate sodium silicate gelation with permeability reduction using the Kozeny-Carman equation (Hiorth et al. 2016(Hiorth et al. , 2017)).However, Eclipse's limitation to communicating only through Satnum files (relative permeability and capillary pressure curves) while running, restricts the modeling of geochemical processes.Therefore, IORSim is currently coupled with SLB Intersect, which allows external programs to update up to 282 static and dynamic simulation parameters.Moreover, coupling IORSim with an open-source reservoir simulator such as OPM Flow would allow glue codes to link the flow equation and geochemical species balance more efficiently together (Rasmussen et al. 2021;Bungartz et al. 2016;Chourdakis et al. 2021).Using reservoir simulators that offer a broader communication interface will allow IORSim to update, for instance, absolute permeability/porosity changes or the tracking and updating of partitioning components in the brine, oil, and gas phases.
The mathematical and numerical concepts of IORSim are introduced in Sect.2, while the IORSim concept is analytically and numerically validated in Sect.3. Section 4 showcases a field-scale reservoir study to demonstrate the seamless integration of IOR-Sim into existing reservoir models (Ekofisk field), the ability to handle comprehensive injection/production histories, and the capability to compute and visualize geochemical reactions in complex geological systems.While the analyzed produced water chemistry shows a clear trend of magnesite and anhydrite precipitation, the presented IORSim study does not aim to make general assumptions about the Ekofisk reservoir.

IORSim Program Structure
The IORSim software provides two dedicated simulation modes that essentially differ in the implemented host simulator and IORSim coupling mechanism.During an IORSim forward mode simulation, the host simulator and geochemical reaction tool coupling is unidirectional.After completing the flow calculation in the host simulator, IORSim is used as a post-processor tool that advects chemical species along the flow direction (Fig. 1a).While the geochemical reactions impact the water chemistry, the global fluid flow remains unaffected.The intentional exclusion of a host simulator feedback mechanism allows high-speed chemical species balance calculations in complex reservoir systems.
If mineral dissolution, precipitation, or other geochemical mechanisms affect the reservoir fluid flow, IORSim provides a backward coupling mechanism.Unlike the IOR-Sim forward mode, the backward coupling mode communicates with the flow simulator at each time step.A sketch of the bidirectional IORSim backward mode is depicted in Fig. 1b.After initiating the reservoir model and calculating the first time step in Eclipse, a Python control script pauses the reservoir flow simulation software (Fig. 1b, II).Since the flow calculations are only paused and a re-initialization of the model is avoided, the additional computation costs are kept to a minimum.IORSim reads the cell phase flow volumes (printed to root.unrst) and well rates (printed to root.rft) to detect the global fluid flow and to permute the grid cells into a sequence of upstream and downstream cells.Instead of assembling and handling a global matrix, the topological sorting algorithm allows IORSim to solve the reactive transport species balance as an iterative string of one-dimensional problems.Section 2.1 and Appendix 6 present a detailed description of the implemented sorting algorithm.After solving the geochemical reactions, IOR-Sim calculates relative permeability and, optionally, capillary pressure curves that are communicated back to Eclipse by using the Eclipse Satnum option (Fig. 1b, III).The continuous process of starting and pausing the Eclipse and IORSim simulation is continued until the reactive transport simulation is completed.Once the last simulation step is reached, the Python control script merges the Eclipse and IORSim results into a single root.unrstfile, which allows standard post-processing application programs such as OPM ResInsight (Dale et al. 2012) or SLB FloViz to visualize and analyze the IORSim reactive transport simulation results.
The IORSim computation performance significantly varies based on the selected (forward or backward) simulation mode and the selected geochemical module.For instance, the Ekofisk IORSim forward simulation presented in Sect.4.2 required around 1 h of computation time on an i9-12th generation processor with 811 time steps, Δt = 20 day, and 56,498 active cells.On the other hand, an IORSim backward mode water diversion study (simpler geochemical module) on the NCS Snorre field required approximately 10 h with 2043 time steps, Δt = 1 day, and 512,550 active cells.

Topological Sorting: From 3D to 1D
Regardless of whether the IORSim forward mode or backward mode is selected, the software uses a sorting algorithm to simplify the transport and reaction calculations of geochemical components.After reading the flow field from the host simulator, IORSim permutes the grid cells with respect to the fluid flow direction.The implemented sorting algorithm leads to a generalized streamline approach in which the original three-dimensional problem is transformed into a sequence of smaller one-dimensional problems.Instead of assembling a global matrix to solve for all unknowns in all grid blocks at once (Fig. 2, I), IORSim detects the fluid flow and updates the geochemistry in one grid block at a time.IORSim initially calculates the geochemical balance at the injection well and then continues to solve the species balances along the fluid flow (Fig. 2, II).Since the species concentration in the depicted center block depends on three neighboring cells, IOR-Sim takes a "detour" to compute the chemical balances in cells 2, 4, and 6 before solving the chemical species balance in cell 7 (Fig. 2, III).After resolving the inflow terms of the center block ( ∑ f =2,4,6 F in f (t + Δt)Δt ), IORSim calculates the outflow term ( F out 7 (t + Δt)Δt ) and continues to solve the chemical species until the production well is reached.
Formally, the water flow field is viewed as a directed graph (Diestel 2016): A directed edge connects block i to block j whenever water flows from i to j.In the absence of flow cycles in the flow field, the graph is a directed acyclic graph (DAG).In this case, a simple depth first search (DFS) algorithm is used to order to the reversed flow graph along the flow direction.
The geochemical component concentrations at the next time level C n+1 are calculated by using Eq. 1, where V w = S w ΔxΔyΔz is the volume of water in the grid cell, C n is the old concentration of the given component, Δt is the time step, F in is the incoming flow of water, Q in is the incoming flow of geochemical components, and r is a source term to account for mineral precipitation/dissolution. The incoming water F in and geochemical component flow Q in is defined by the sum over the neighboring grid blocks and positive by definition A derivation of the chemical species balance (Eq. 1) implemented in IORSim is provided in Appendix 6.
The authors acknowledge that in addition to the cell permutation technique used in IORSim, Natvig, Lie, and colleagues have independently developed a topological sorting method (Natvig and Lie 2008;Lie et al. 2014).Moreover, Natvig and Lie (2008) and Lie et al. (2014) discussed the possibilities of generalizing the sorting approach to cases in which the flow field contains closed flow cycles.After identifying strongly connected components (SCC) inside a flow graph, Natvig and Lie suggested solving each subsystem with a nonlinear solver before updating dependent grid cells.Several well-known algorithms exist for finding the SCC of a directed graph, e.g., Tarjan's algorithm (Tarjan 1972) or Kosaraju's algorithm (Sharir 1981). (1)

1D Tracer Advection
In addition to a fully implicit species transport scheme, IORSim provides an explicit and adaptive numerical scheme.The selection of the numerical scheme is governed by the Courant-Friedrichs-Lewy (CFL) condition, in which a zero CFL criterion leads to a purely implicit solution and a value of one results in a strictly explicit solution (Courant et al. 1928).
The chemical species transport implementation was validated against an analytical transport calculations of a one-dimensional tracer core flood as illustrated in Fig. 3a.After four days of pure water injection, a step-pulse of a non-partitioning tracer is injected.In general, tracer transport is governed by the advection-diffusion equation with constant diffusion D where the first term of the right-hand side denotes the advective flow, the second term denotes the chemical reactions, and the third term denotes mechanical dispersion.While the injected tracer concentration front follows the average water flow velocity through the column, the heterogeneous pore network causes variations in the flow velocity.Instead of obtaining piston flow, the non-uniform flow velocity causes a (mechanical) dispersion of the tracer concentration front.Assuming a semi-infinite column and a constant injection rate, the effluent tracer concentration c f can be calculated as follows (Appelo and Postma 2004;Lapidus and Amundson 1952): where c inj is the injected tracer concentration, L is the total core length, v is the flow veloc- ity, D L is the longitudinal dispersion, and erfc is the complementary error function.Cur- rently, IORSim does not consider mechanical dispersion, but the numerical scheme introduces numerical dispersion.Assuming constant Δx and Δt , Lantz (1971) derived numerical dispersion coefficients for the explicit and implicit schemes where x is the cell length, Δt is the time step size, and v is the flow velocity.
The implicit and explicit analytical and IORSim tracer transport solutions for a time step size of 0.1 days are plotted in Fig. 3b, c.For t = 0.1 days, the cell inflow at one time step is equivalent to the cell volume, which corresponds to a CFL criteria of one.Figure 3b shows that the explicit numerical dispersion term becomes zero, and the tracer breakthrough concentration approaches one.The zero dispersion case represents a hypothetical case in which the core material consists of a perfectly uniform porous media.While the explicit solution shows a piston-like species movement, Fig. 3c demonstrates that the implicit tracer front is significantly smeared out.Shorter time steps are required to decrease the impact of numerical dispersion on an implicit transport scheme.Figure 3d, e displays the analytical and IORSim tracer transport solution for a time step size of 0.01 days.In this case, the implicit and explicit tracer transport solutions are characterized by a similar (small) numerical dispersion.
The analytical tracer transport calculations match and validate the implicit and explicit IORSim schemes.Moreover, the numerically calculated concentration/time integrals align with the mass injected in all four plotted tracer experiments, confirming mass conservation.

3D Tracer Advection
Section 3.1 analytically validates the implicit and explicit IORSim chemical species transport.However, while one-dimensional tracer transport represents the most straightforward flow possible, geological formations are typically represented by unstructured grids.Most finite difference reservoir simulators use a seven-point stencil scheme to solve the partial differential equation of a cell and its six neighboring cells (Ron et al. 2007).Moreover, Eclipse 100 and other reservoir simulators support explicitly formulated cell connections to model diagonal flow along faults or pinched-out regions.Therefore, in a realistic reservoir case, the IORSim sorting algorithm must correctly identify the flow of a cell and its six surrounding cells, consider possibly pinched-out grid blocks, and integrate non-neighboring cell connections. (5) A reservoir segment of the NCS field Snorre was selected to test the IORSim chemical species transport solution on a geologically realistic field case.Besides diagonal flow along faults, the Snorre reservoir model contains several pinched-out reservoir zones.In the implemented species transport comparison example, 125 kg of non-partitioning tracer mass was admixed to the injection fluid.The injected tracer predominantly migrated along a thief-zone and was produced through a single production well.
The blue points in Fig. 4a represent the tracer production concentration observed using a standard Eclipse 100 simulation.The black curve displays the non-partitioning tracer production signal using a classical Newton-Raphson solver in IORSim.Finally, the red line illustrates the tracer transport profile with the enabled IORSim sorting algorithm.

Numerical Dispersion
The literature discusses several approaches to reduce numerical dispersion.Firstly, Lantz (1971) showed that explicit schemes are less impacted by numerical dispersion than implicit schemes.While numerical dispersion is proportional to x − Δtv in an explicit scheme, truncation errors in implicit schemes are proportional to x + Δtv (cf.Eqs. 6and Fig. 4 The three-dimensional IORSim local implicit tracer advection agrees with the IORSim global implicit and Eclipse 100 tracer transport result (a).The performance advantages of the local implicit IOR-Sim solver increases with increasing number of components (c).Compared to the selection of an implicit and explicit species transport computation, grid refinement causes the strongest numerical dispersion reduction (b).Using splay trees to retrieve chemical equilibriums can reduce computation times by up a factor of up to ten, while yet providing accurate numerical solutions (d) 7).Consequently, as shown in Sect.3.1, numerical dispersion in implicit schemes increases with increasing time steps.While high-order discretization schemes are acknowledged to be significantly less impacted by numerical dispersion, high-order discretization schemes require higher computation efforts than upwind finite volume discretization (Mykkeltvedt et al. 2017).
As the species transport calculations are performed on a duplicated (IORSim) grid, IORSim provides a highly efficient and fast grid refinement option to limit chemical species dispersion.For a regular non-well reservoir grid block, the refined velocities are calculated by spatially weighted interpolation between the Eclipse flow rates on each side of the refined block connection.The IORSim grid refinement options allow the definition of one or several finer grid resolution regions with equal or different degrees of refinement.Figure 4b compares the produced tracer concentration when selecting an implicit, explicit, and refined IORSim chemical species transport scheme.In line with the literature, the implicit solution shows the strongest numerical dispersion when testing non-reactive tracer transport on the Snorre grid.The most substantial decrease in numerical dispersion was observed when enabling the IORSim grid refinement.

Sorting Algorithm Efficiency
The IORSim sorting algorithm was developed to address two main challenges of reactive transport simulations.Firstly, ensuring numerical stability while solving nonlinear geochemical reactions and, secondly, reducing computation times.
At an early stage of the software development, the IORSim tracer module used a classical, fully implicit approach to solve the chemical species balance.To demonstrate the computation efficiency of the IORSim sorting algorithm, the computation times of a global and a local implicit tracer simulation scenario were compared.Figure 4c shows the normalized IORSim computation times as a function of increasing non-partitioning tracer components.While the global and local implicit solver perform similarly when a small number of tracer species is injected, the efficiency of the IORSim sorting algorithm becomes distinct with an increasing number of tracer species.Compared to the global Newton-Raphson solver, the transport simulation of 100 tracer species (176,664 active cells) was approximately four times faster when activating the IORSim sorting algorithm.Although the non-reactive transport computation comparison certifies IORSim's high-performance efficiency, the IORSim sorting algorithm is expected to become more significant when calculating nonlinear geochemical reactions.

Accelerating Geochemical Computations
IORSim offers a wide range of geochemical modules with varying levels of complexity.As the non-partitioning tracer, partitioning tracer, and sodium silicate gelation modules solve a limited set of nonlinear reactions, the decoupled chemical species balance approach is sufficient to obtain reasonable computation times.However, if IORSim addresses more intricate modeling tasks, such as the visualization of the initial and injected brine composition impact on reservoir chemistry (Sect.4.4), IORSim must solve interconnected geochemical reactions (such as mineralization, ion exchange, surface complexation, etc.) of multicomponent systems.In this case, IORSim applies basis swapping to reduce the number of unknowns and utilizes splay tree calculations to reuse geochemical equilibriums that have been solved before.

Basis Swapping
The IORSim solver divides the geochemical species into a set of basis and secondary species to speed up the calculation (Morel 1983;Bethke 2022).Based on the vector spaces theory in linear algebra, the concentrations of secondary species are uniquely determined from the basis species concentrations.As a result, IORSim reduces the number of unknowns by only solving for the basis species concentrations directly.While, in principle, any linearly independent set of species of the right dimension is suitable to determine the species compositions, IORSim defines the major reservoir ions and proton concentrations (i.e., the "free ion form" of each element) as the basis species.Moreover, in cases where a solution is assumed to be in equilibrium with one or more minerals, IORSim removes one of the major ions that the mineral consists of and replaces it with the mineral.The implemented basis swapping technique is essentially the same as presented by Bethke (2022).

Reusing Old Solutions with Splay Trees
While the IORSim sorting algorithm and basis swapping reduce computation time, complex geochemical simulations are yet computationally demanding.The literature describes the application of proxy models, e.g., various purely statistical algorithms, and/or machine learning approaches to reduce computation times (Jatnieks et al. 2016;Demirer et al. 2022; Laloy and Jacques 2022).However, as machine learning techniques tend to introduce mass balance errors (Guérillot and Bruyelle 2020), other research groups suggested solving chemical equations while reducing the precision of the solution.For example, Leal et al. (2020) and Kyas et al. (2022) combined an on-demand extrapolation strategy based on first-order Taylor approximations to speed up the computations by one to three orders of magnitude.De Lucia et al. ( 2021) focused on parallelization and the incorporation of distributed hash maps to store and retrieve solutions of previously solved geochemical problems.While similar-looking geochemical problems are treated equally, the interpolation or extrapolation of geochemical reactions is avoided.
IORSim uses a splay tree method to reduce the number of geochemical reaction calculations (Hiorth et al. 2016(Hiorth et al. , 2017)).If the geochemical solver encounters an equilibrium or an identical/similar problem solved before, IORSim retrieves the geochemical solution from splay trees (Sleator and Tarjan 1985).Initially, the program stores the geochemical solution concentrations in a splay tree based on the input concentration, with a new splay tree being initialized for each combination of minerals.If the grid cell concentrations are close enough to a previously calculated solution, the new species solution is obtained from the splay tree.A unique aspect of the IORSim implementation is the creation of splay tree keys and the retrieval of stored data.While only aqueous species concentrations are used to construct lookup keys, equilibrium calculations determine the adsorbed concentrations.Figure 4d demonstrates that the implemented IORSim splay tree technique can reduce computation times by a factor of up to ten while providing accurate numerical solutions.Compared to the computation time of the reference case (in which the chemical equilibriums are continuously solved over iterations), the computation time reduces with decreasing splay tree resolution.A detailed explanation of the implemented splay tree lookup algorithm is provided in Appendix 7.

Geochemical Solver Validation
The IORSim geochemical solver incorporates the SUPCRT92 database of Johnson et al. (1992) to calculate the thermodynamic properties of minerals, gases, aqueous species, and reactions along characteristic subsurface temperatures and pressure conditions.Details on the implemented aqueous equilibrium chemistry, surface chemistry, and mineral kinetics are provided in Appendix 8.1 to Appendix 8.2 Before implementing a full-scale reservoir Ekofisk study, the geochemical IORSim solver was validated against (a) the well-documented PHREEQC ion transport and exchange benchmark (Parkhurst and Appelo 2013;Moortgat et al. 2020;Kolditz et al. 2012;Rolle et al. 2018) and (b) a core flooding experiment reported by Hiorth et al. (2013).
The PHREEQC benchmark represents a column saturated with a sodium-potassiumnitrate solution in equilibrium with a cation exchanger.Once a calcium chloride solution is injected into the sample, the calcium, potassium, and sodium ions react with the exchanger until chemical equilibrium is reached.Since chlorine acts as an inert ion, the chlorine effluents solely depend on the transport term.Figure 5a compares the IORSim geochemical solver with the PHREEQC solver results.While the overall effluent production is excellent, minor concentration differences are explained by the differences in underlying geochemical databases and the differences in employed numerical schemes.
During the second geochemical solver benchmark test, the IORSim module was tested against experimental core flooding data reported by Hiorth et al. (2013).The selected Liège outcrop sample was initially saturated with deionized water and then flooded by a calcium, magnesium, and sulfate-rich brine (Table 2).Figure 5b depicts the ion effluent concentration observed in the single-phase flooding experiments conducted at a temperature of 130 • C.After two days of injection, the effluent chlorine concentration matched the injected concentration, indicating a complete displacement of the deionized water (Table 2 and Fig. 5b, yellow curve).While the inert chlorine ions do not participate in geochemical reactions (Appelo and Postma 2004), the concentration differences between the injected and produced divalent ions reveal geochemical reactions.Three geochemical processes were considered to reproduce the divalent ion concentration profile: ion exchange, ion adsorption, and mineral kinetics.
Ion exchange is a rapid process that involves ion swapping between the solution and the mineral surface without altering the mineral's structure.In Fig. 5b, ion exchange causes the calcium concentration to peak approximately two days after the brine injection starts.In contrast to ion exchange, the more gradual ion adsorption involves ion adhesion from a solution onto the pore surface without necessarily releasing other ions back into the solution, hence potentially changing the surface charge.The effluent ion concentration after ten simulation days (Fig. 5b) is predominantly impacted by the simulation's mineral kinetic rates.Mineral dissolution/precipitation is a slow geochemical process that involves the breakdown or formation of the mineral's crystal structure.Nermoen et al. (2015) experimentally demonstrated that magnesite precipitation in chalk samples can occur over several years.
The numerical IORSim effluent results are depicted in Figure 5b and the utilized rate mineral precipitation/dissolution constants, ion exchange rate, and surface complexation equilibrium rates are listed in Table 1.While IORSim replicates the effluent trend, constant mineral kinetics rates cannot reproduce the experimentally obtained slow-reaching mineral equilibrium.A comprehensive discussion and suggestion of dynamic mineral equilibrium constants are provided in Pedersen et al. (2016).

Ekofisk Field
Under production since 1971, the Norwegian offshore field Ekofisk is considered the first major oil discovery in the North Sea.Separated by an impermeable layer, an upper (Ekofisk) and lower (Tor) formation comprise the chalk reservoir.While the Ekofisk formation is highly fractured, the matrix porosity varies between 25 and 40 %, and the permeability ranges from 0.1 to 5 mD.Especially in the upper Ekofisk formation, mineralization, unfractured chalks, and tight clays act as local flow barrier (Jakobsson and Christian 1994).
The Ekofisk anticline structure spans an area of approximately five by ten kilometers and can be divided into smaller segments to allow faster and more efficient reservoir modeling.Figure 6 2).
Although the Ekofisk formation contains silica, illite, smectite, and kaolinite tracers, the upper Ekofisk and the lower Tor structure are considered pure calcite reservoirs (Kallesten et al. 2021).Numerous studies have investigated the brine-rock interactions caused by seawater injection into chalk formations.Madland et al. (2011) demonstrated that seawater injection into Liège and Stevns Klint chalk samples resulted in the precipitation of anhydrite and magnesium-bearing minerals.Similarly, Megawati et al. (2015) showed that the injection of magnesium-rich brines caused calcite dissolution and promoted magnesite precipitations.X-ray diffraction analysis (XRD) confirmed that magnesite precipitations occurred in pure as well as impure chalk samples.Moreover, Minde et al. (2018) highlighted the impact of temperature on calcite dissolution and magnesite precipitation.If magnesium-rich brines were injected into chalk samples, energy-dispersive X-ray (EDX) and XRD measurements detected magnesite crystals in flooding experiments that exceeded 110 • C.

Non-Reactive Transport Ekofisk Simulation
To visualize and quantify the geochemical interactions induced by seawater injection on the Ekofisk field, an IORSim forward mode simulation was implemented.The field operator supplied the Ekofisk Eclipse 100 reservoir model (comprising 89,936 grid blocks), covering the southern segment of the Ekofisk formation.While the reservoir model included 53 wells, the IORSim modeling study focused on (1) replicating the historical P-1 production rates, and (2) matching the corresponding produced water composition.Figure 6 marks the P-1 location in the lower left corner of the reservoir model, which communicates with the seawater injection well I-1.Separated by approximately 400 ms, the seawater injection well I-1 penetrates both the upper Ekofisk and lower Tor formations, while the production well P-1 is perforated in the upper Ekofisk formation.The combined effects of pressure depletion and rock-fluid interactions have led to significant Ekofisk reservoir compaction, which is numerically captured by the Eclipse 100 rock compaction module.To benefit from the history-matched Ekofisk model provided by the field operator, the presented IORSim Ekofisk study was run in forward mode, meaning that the chemical components were advected along the reservoir flow after the Eclipse simulation was completed.However, to further improve the accuracy of field history-matching, it would be feasible to conduct an IORSim backward mode study that correlates mineral dissolution and precipitation with changes in permeability and porosity.
After reviewing the field history matched, the sampled chlorine ion concentrations were used to determine the formation water and seawater mixing proportion of well P-1.Figures 7 and 8a show the water volume and chlorine concentration history-matches that were obtained by (1) balancing the fracture permeability introduced between injector I-1 and producer P-1, and (2) fine-tuning the well transmissibilities of I-1 and P-1.
As an inert ion, the chlorine concentration depends solely on the formation and injection water mixing.Between 1997 and 2000, the produced water composition agreed with the Ekofisk formation water.A sharp chlorine concentration decline from 2000 onwards indicates the seawater breakthrough at the production well.After two years of decreasing chlorine concentration, both sampled and modeled chlorine levels aligned with seawater chlorine concentration (Fig. 8a).
In contrast to chlorine, the concentration of the reactive calcium, magnesium, and sulfate ions deviated significantly from the measured data if geochemical reactions were not considered.Figure 8b and Table 2 show that ten years after seawater breakthrough (2010), the produced P-1 water calcium concentration was approximately 0.02 mol/L higher than the calcium concentration of the injected seawater.At the same time, the produced water chemistry is characterized by an approximately 0.03 mol/L lower magnesium and 0.02 mol/L lower sulfate ion concentration than the injected seawater (Fig. 8 and Table 2).

Rate-Controlled Mineralization
The significant mismatch of the non-reactive divalent ion transport and sampled ion concentration confirms that seawater injection induces chemical reactions.Activating mineral kinetic reactions in IORSim led to an improved match of the simulated and measured effluent ion concentration, as demonstrated by the green curves in Fig. 8b-d.The chemistry of seawater and brines is very different from the chemistry of rivers and drinking water; for a comprehensive discussion on this topic, see Chapter 4 in Garrels and Christ (1965).In Appendix 8.4, we have listed the chemical reactions included in IORSim used to describe seawater chemistry.When seawater, which at ambient conditions is close to being in equilibrium with calcite, is heated up to reservoir temperatures ( ∼ 130 • C), magnesite becomes supersaturated.Magnesite will precipitate, and since Fig. 8 Produced water composition from well P-1.After using inert ion concentration to replicate the seawater breakthrough, the non-reactive calcium, magnesium, and sulfate ion transport simulations show a substantial mismatch between the simulated and measured divalent ion concentrations (blue curves).When the rate-controlled mineralization (green curves) and surface complexation (red curves) modules were activated in IORSim (reactive transport simulation), the produced water chemistry match was significantly improved magnesite contains carbonate, carbonate is removed from the solution.When carbonate is removed, calcite becomes under-saturated and dissolves.The net effect of these reactions is summarized in the equation as follows: Thus, magnesite precipitation leads to an increase in the calcium concentration in the solution.When the calcium concentration increases, anhydrite (calcium sulfate) becomes supersaturated and will precipitate.
The mineral dissolution and precipitation kinetics in IORSim are implemented according to Palandri and Kharaka (2004) and are presented in Appendix 8.2.

Surface Complexation
While the rate-controlled mineral kinetics capture the trend of the reactive ions concentration changes, particularly the sulfate ions breakthrough delay indicates the presence of additional geochemical reactions.Besides mineral kinetics, the second reactive transport IORSim Ekofisk simulation used surface complexation to improve the produced water chemistry match (Fig. 8b-d, red line).The sulfate adsorption is implemented as follows: The IORSim surface complexation model was developed by Van Cappellen et al. (1993) and Hiorth et al. (2010) and is summarized in Appendix 8.3.

Conclusion
IORSim is a reactive transport simulator designed as a plug-in tool for industry-standard reservoir flow simulators.While the host simulator is used to integrate reservoir models and to perform the global flow calculations, IORSim solves the decoupled chemical species balance.The novel software coupling approach allows the implementation of fullfield reservoir simulations with a large number of grid cells and chemical components.
• The IORSim sorting algorithm decouples the nonlinear geochemical reaction calculations into a sequence of recurring one-dimensional problems to assure numerical stability and computation efficiency.( After introducing the Eclipse/IORSim coupling mechanism and deriving the mathematical IORSim concept, the implemented program structure was analytically and numerically verified on a set of demonstration cases.IORSim was then used to investigate the impact of seawater injection on water chemistry and mineralogy in the NCS field Ekofisk.
• The decoupled flow simulator and geochemical reaction calculations allow seamless integration of full-field reservoir models that contain complex geological structures, a large number of wells, and long production histories.• The reactive IORSim modeling of the Ekofisk case demonstrates that seawater injection into chalk formations impacts the mineralogy and changes the produced water chemistry.In the investigated Ekofisk case, seawater promoted calcite dissolution, and led to the precipitation of magnesite and anhydrite.Moreover, surface complexation modeling revealed that sulfate is adsorbed on the calcite surface.

Appendix 1: Sorting Algorithm Derivation
IORSim uses a cell-centered finite volume method (LeVeque et al. 2002).For a fixed volume Ω , the integral form of the mass balance equation for an arbitrary component is where Γ = Ω denotes the volume boundary, is the reservoir porosity, C is the molar concentration of the chemical component, n is the outward-pointing normal vector to the boundary, J is the total mole flux of the component, and r is the source term.IORSim uses an implicit scheme to integrate Eq. 11 forward from t n to t n+1 = t n + Δt .For simplicity, the derivation of Eq. 1 is shown for the case of one-dimensional flow in a Cartesian grid with a constant grid size, Δx , and a constant cross-sectional area, A. However, the derivation can be extended to the general cases.By abuse of notation, the concentration of the chemical component C and the source term r are initially replaced by their volume-averages.The total volume of each grid block is then Ω = AΔx , and it is straightforward to derive the discretized Eq. 12 where C i−1∕2 is the concentration at the left boundary of the grid block, C i+1∕2 is the con- centration at the right boundary, and Q w = uA denotes the volumetric flow rate of water across the interface between two cells.
So far, this is a conventional implicit scheme that can be solved as sparse system of coupled, nonlinear equations.The critical point is to rewrite the system so that the concentration in block i at time t n+1 solely depends on (1) the concentration in block i at the previous time and (2) the concentrations in the upstream cells that are already updated to the next time level.To this end, it can be assumed without loss of generality that Q w > 0 .In this case, the concentration at the left boundary is also the inflowing concentration, C in = C i−1∕2 = C i−1 .Similarly, the concentration at the right boundary then becomes the outflowing concentration, C out = C i+1∕2 = C i .Equation 12 can therefore be written as where and In the above equations, V w = S w AΔx is the volume of water.For incompressible flow, Q w,in = Q w,out , which yields where B ≡ Q w Δt∕V w was introduced.The physical interpretation of B is that it is the frac- tion of the grid block being flooded during the time step.When applying the exact same procedure to the case of variable grid block dimensions and a three-dimensional grid, the derived equations are the same, provided each of the two flux terms is replaced with a sum to account for the possibility of multiple neighbors The first sum is taken across all cell interfaces with an incoming flux ( Q in,f > 0 ), and the second sum over all interfaces with an outgoing flux ( Q out,f > 0 ).Thus, by definition, all flow rate terms remain positive in the notation.Hiorth et al. (2017) noted the possibility of small mass balance error in the host simulator that impact the IORSim calculations.Therefore, the finite volume formula was corrected to take the volume balance of the flow simulator into account where In other words, the denominators of Eqs. 14 and 15 are replaced by inserting Eq. 19, which leads to Eq. 1 listed in Sect. 2.

Appendix 2: Splay Tree Lookup Algorithm
To save time in complex full-field simulations, IORSim uses splay trees to reuse solutions to geochemical equilibrium problems.However, due to the appearance or disappearance of minerals, it is more complex to store and search for solutions.IORSim creates one splay tree object for each possible combination of minerals.Whenever the geochemical solver is invoked, IORSim detects the minerals present at the given point in space and time.The information is encoded in a binary key, which is used to recover the correct splay tree.For each splay tree, the following pieces of data are used to construct lookup keys: • The total concentrations of each basis species (including surface basis species).
• The temperature, T.
• The log activity of minerals.
Given this data, the equilibrium state is uniquely determined.Each node in a splay tree stores the following pieces of data: • Differences in the total concentration of each basis specie over the time step (old values minus new values).• If ion exchange and/or surface complexation reactions are included, differences in the adsorbed concentrations (old-new).• If surface complexation reactions are included, (surface charge density) and (surface potential).• The pH value.
To be more specific, each key consists of n B + n M + 2 integers, where n B is the maximal number of basis species that can exist in the entire simulation run, and n M is the maximal number of minerals.The first n B + 1 entries are computed from the basis specie concen- trations and the time step: given an input value v, it is converted into an integer via the formula where N int is the input splay tree resolution, f log is a scaling factor (here, assumed to be unity), and the brackets signifies a floor-function.The splay tree resolution is related to w via On the other hand, the final n M + 1 values are obtained from the temperature and the min- eral log activities.Each input value v is converted in to an integer key using Eq.23 where f lin is a scaling factor that is set to unity.Essentially, the approach categorizes solu- tions into bins; the finer the resolution of the splay tree, the greater the number of bins ( (meaning, more solutions are considered distinctly).When a new scenario arises, the solver initially seeks the "closest" pre-computed solution by utilizing the key.IORSim enables the user to adjust the splay resolution by setting an interpolation parameter.While a finer splay resolution causes longer computation, numerical accuracy increases with larger interpolation parameters.
Equations 27-33 are solved together with the reactions in the aqueous phase and are thus included in the mass balance equations.In the absence of data to the contrary, it is assumed that the temperature-dependence of the surface reactions is the same as for the analogous reactions in the aqueous phase (Hiorth et al. 2010).
A key difference between the aqueous and surface reactions is that the activity at the surface is related to the surface potential, 0 [V].If the surface is negatively charged, positive ions will have a higher activity compared to negative ions, while the opposite is valid for a positively charged surface.
The relationship between the surface potential and the surface charge density, [ C∕m 2 ], is given by the Grahame equation Equation 34 is also included in the set of nonlinear equations to solve in each grid block at each time step; the physical solution is where and 0 have the same sign.
• denotes the relative permittivity (dielectric constant) of water.• R = 8.314 J/K/mol is the ideal gas constant.
• T is the absolute temperature.
• n B i is the molar concentration of charged species number i in the bulk aqueous phase away from the mineral surface.
Finally, the surface charge density is obtained from where S A is the specific surface area in units of square meter per liter of water, and the brackets denote molar concentrations of the surface complexes.For details on the numerical implementation, see Nødland and Hiorth (2022).(31) > CO 3 H 0 ⇌ > CO − 3 + H + .

Fig. 1
Fig. 1 IORSim and host simulator coupling schemes.The IORSim forward mode utilizes unidirectional coupling to allow an efficient and fast visualization of chemical species and mineralogy distribution (a).The impact of geochemical reactions on global fluid flow is calculated when selecting the bidirectional IORSim backward mode (b)

Fig. 2
Fig.2Chemical species balance.After sorting the cells along flow lines from the injection well to the production well (sorted cell numbering (II)), the chemical species balance is solved as a sequence of (nonlinear) one-dimensional problems (III)

Fig. 3
Fig. 3 Analytical IORSim validation.The implicit and explicit IORSim tracer concentration profiles show an excellent agreement with the corresponding analytical solutions.The analytical and numerical solution assume a core length L of 0.3048 m, flow velocity of 0.3048 m/day, and cell length X of 0.03048 m

Fig. 5
Fig. 5 IORSim reproduces the PHREEQC ion transport and case (Parkhurst and Appelo 2013).The IORSim geochemical solver replicates the numerical and experimental brine effluent composition results reported by Hiorth et al. (2013) (b) displays the southern Ekofisk reservoir model in which 19 injection and 34 production wells were operated between Mai 1975 and June 2014.The reservoir model includes three distinct brine compositions.Besides the high saline Tor formation brine (TDS of 192.7 g/L) and the medium saline Ekofisk formation brine (TDS of 83.5 g/L), around 50 Mm 3 seawater was injected into the reservoir (TDS of 36.1 g/L) (Table

Fig. 6
Fig. 6 Southern Ekofisk reservoir model with wells operated in April 1997.The production well P-1 is located on the lower left part of the reservoir model and communicates with the seawater injection well I-1

Fig. 7
Fig. 7 Historical and modeled P-1 produced water rates • Moreover, for a complex geochemical problem, IORSim provides basis swapping and splay trees techniques to reduce computation times significantly.• IORSim offers two dedicated simulation modes: The high-speed IORSim forward mode allows the post-advection of chemical components to visualize species distribution, water chemistry, and mineral interactions.If the geochemical reactions interfere with the fluid flow, the IORSim backward mode uses relative permeability to update the global fluid flow at each time step.(8)Mg 2+ + CaCO 3 (s) ⟶ Ca 2+ + MgCO 3 (s).

Table 1
Parameters used to reproduce the core flooding ion effluent concentration depicted in Fig.5bThe mineral kinetic constants are denoted by r k and E a is the activation energy.Equilibrium constants are listed in Appendix 8.4

Table 4
Table of computed log K values for surface and aqueous complexes in the carbonate surface complexation modeling