A Transient Boundary Element Framework for Pressure Buildup, Injectivity Enhancement, and Reactive Permeability Evolution During CO₂ Storage in Heterogeneous Carbonate Reservoirs

2.1. Mass conservation of the pore fluid

The fluid mass balance in a deformable porous medium may be written in the form

To account for the compressibility of both fluid and porous skeleton, it is convenient to introduce a generalized storage coefficient $S$ , so that the pressure equation can be written as

For a point or line injection source, the well may be represented by an internal singular source term,

2.2. Darcy flow law with evolving hydraulic conductivity

The Darcy velocity is governed by

where $k(\mathbf{x},t)$ is the intrinsic permeability, $\mu$ is the dynamic viscosity, and $\mathbf{g}$ is the gravitational acceleration vector. For horizontal injection problems or when pressure perturbations dominate over gravity, the body force term may be neglected, yielding

Substituting Darcy’s law into the mass balance equation gives the governing hydraulic diffusion equation

This equation is nonlinear because both $\varphi$ and $k$ evolve with time as a result of reactive dissolution and stress redistribution. In a heterogeneous reservoir, the permeability field may also vary spatially from one subregion to another, such that

2.3. Reaction kinetics for carbonate dissolution

The geochemical alteration of the reservoir is represented through a scalar reaction progress variable $\xi (\mathbf{x},t)$ , where $\xi =0$ corresponds to the undegraded material state and increasing $\xi$ indicates progressive mineral dissolution. At the Darcy scale, the kinetics may be expressed in a reduced phenomenological form as

This form captures three physically relevant features: reaction is initiated by injection-induced disequilibrium, the rate increases with pressure perturbation, and the process gradually saturates as $\xi \to 1$ . It also provides a convenient scalar internal variable for coupling geochemical alteration to porosity and permeability evolution.

More elaborate kinetic laws can also be used without changing the overall mathematical framework. For example, $\mathrm{\Psi }\left(p\right)$ may be replaced by a function of dissolved CO₂ concentration or mineral saturation. In this study, the pressure-controlled form is adopted because it captures the essential reactive behavior while keeping the formulation efficient and compatible with the boundary element method.

2.4. Porosity evolution law

The dissolution of carbonate minerals enlarges pore space and changes the local storage characteristics of the medium. This effect is incorporated through a constitutive relation linking porosity to the reaction progress variable and, if desired, to volumetric strain. A general differential form is

For the present carbonate dissolution model, the dominant chemical contribution is represented as

Substituting the reaction law into the porosity relation yields

2.5. Permeability evolution law with reactive and stress-sensitive effects

Permeability evolution is one of the most important nonlinear mechanisms in reactive CO₂ injection because it controls both injectivity and the spatial redistribution of pressure. Experimental observations in dissolving carbonate formations often indicate a highly nonlinear increase of permeability with pore enlargement. To capture this effect, the present model adopts a generalized constitutive relation of the form

A convenient and physically meaningful representation is the multiplicative law

If a purely chemical model is desired, the permeability law reduces to

Alternatively, a Kozeny–Carman-type relation may be used:

2.6. Geomechanical equilibrium and effective stress coupling

To incorporate geomechanical feedback, the porous solid is modeled as a linear poroelastic continuum under quasi-static conditions. Neglecting inertia, the equilibrium equation is

For small strains, the strain tensor is

The volumetric strain relevant for porosity evolution is then

To couple mechanics back to flow, the mean effective stress is defined as

2.7. Strongly coupled pressure–reaction–deformation system

Combining the preceding relations, the governing reactive–geomechanical flow model may be written as the following coupled system in $\mathrm{\Omega }\times (0,T]$ :

2.7.1. Pressure diffusion equation

2.7.2. Reaction progress equation

2.7.3. Porosity evolution equation

2.7.4. Geomechanical equilibrium equation

2.7.5. Permeability constitutive law

This coupled system is nonlinear through three mechanisms: the dependence of reaction rate on pressure, the dependence of porosity on reaction progress and deformation, and the dependence of permeability on both reaction and effective stress. These nonlinearities are responsible for the feedback processes observed in dissolving carbonate reservoirs, including near-well permeability enhancement, pressure redistribution, and mitigation of long-term pressure buildup.

2.8. Initial and boundary conditions

To complete the mathematical statement of the problem, appropriate initial and boundary conditions must be specified.

2.8.1. Initial conditions

At $t=0$ , the reservoir is assumed to be in a reference equilibrium state:

2.8.2. Hydraulic boundary conditions

The external boundary $\mathrm{\Gamma }$ may be decomposed into prescribed-pressure and prescribed-flux parts, such that

On ${\mathrm{\Gamma }}_p$ , Dirichlet conditions are imposed:

On ${\mathrm{\Gamma }}_q$ , Neumann conditions are prescribed:

For an infinite or semi-infinite reservoir idealization, a far-field condition is adopted:

2.8.3. Mechanical boundary conditions

Similarly, the mechanical boundary is partitioned into displacement and traction boundaries:

On ${\mathrm{\Gamma }}_u$ ,

These conditions allow representation of confinement, overburden loading, symmetry planes, or free boundaries depending on reservoir geometry.

2.9. Interface conditions in heterogeneous reservoirs

Because the reservoir may consist of multiple subregions with distinct initial permeabilities, porosities, and mechanical properties, continuity conditions must be enforced across each internal interface ${\mathrm{\Gamma }}_{\mathit{ij}}$ separating subdomains ${\mathrm{\Omega }}_i\text{and}{\mathrm{\Omega }}_j.$

For hydraulic flow, the continuity of pressure and normal flux requires

For mechanics, continuity of displacement and traction is imposed:

These interface relations are essential for the boundary element formulation because they permit decomposition of a heterogeneous reservoir into piecewise homogeneous or weakly heterogeneous regions while preserving the correct transmission of hydraulic and mechanical fields across internal boundaries.

2.10. Incremental linearization for nonlinear solution

The constitutive relations above render the problem nonlinear, and therefore an incremental solution strategy is adopted. Let $t^{n+1}=t^n+\mathrm{\Delta }t$ . Over each time increment, the variables are updated as

The reaction term is linearized as

This procedure transforms the nonlinear coupled system into a sequence of linearized subproblems that can be efficiently solved within each time step while retaining the evolving hydro-chemo-mechanical feedback. Such an incremental treatment is especially suitable for BEM implementation because it preserves the linear operator structure required for boundary integral transformation.

2.11. Reduced form used in the present study

Although the full formulation above includes explicit geomechanical equilibrium, many practical CO₂ injection problems can be represented by a reduced reactive-flow model in which deformation is absorbed into an effective stress-sensitive permeability law. In that case, the governing equations used in computation reduce to

This reduced system retains the essential physics highlighted in the manuscript, namely transient pressure diffusion, dissolution-driven porosity growth, and nonlinear permeability enhancement near the injection zone, while remaining fully compatible with boundary-only discretization. It therefore forms the foundation of the subsequent boundary integral development and numerical implementation.

3.1. Incremental linearization of the reactive flow operator

At time level $t^{n+1}=t^n+\mathrm{\Delta }t,$ the nonlinear pressure equation in Section 2 is linearized by freezing the coefficients at the previous iteration or previous time level. Let

Over the interval $[t^n,t^{n+1}]$ , the hydraulic operator is approximated using effective coefficients

Introducing the hydraulic diffusivity

This form is the basis for the time-dependent boundary integral representation.

3.2. Fundamental solution for the transient diffusion operator

For each subregion ${\mathrm{\Omega }}_i,$ the corresponding adjoint fundamental solution $G_i(\mathbf{x},t;\mathit{\xi },\tau )$ is defined as the solution of

In two-dimensional infinite space, the transient diffusion fundamental solution is

The normal derivative required in the boundary integral equation is given by

3.3. Time-dependent boundary integral representation

Multiplying the linearized pressure equation in ${\mathrm{\Omega }}_i$ by the fundamental solution $G_i,$ multiplying the adjoint equation for $G_i\mathrm{by}{p}_i,$ subtracting the two expressions and integrating over ${\mathrm{\Omega }}_i\times [t^n,t^{n+1}]$ , one obtains the transient reciprocal identity. After application of Green’s second identity in space and integration by parts in time, the pressure field at a collocation point $\mathit{\xi }\in \overline{{\mathrm{\Omega }}_i}$ satisfies

The preceding expression is exact for the incrementally linearized subproblem. It contains three classes of terms: boundary integrals, internal source terms due to wells and nonlinear corrections, and domain history terms associated with the transient initial state. In classical BEM, domain integrals are absent for homogeneous linear problems; in the present reactive setting, they arise from source forcing and coefficient updates. Their treatment is discussed below.

3.4. Treatment of injection wells and equivalent domain terms

A principal advantage of the present formulation is that physical wells can be represented analytically as internal singular sources without introducing local volumetric refinement. For a well located at ${\mathbf{x}}_w$ , the source term is written as

Thus, the well contribution is evaluated semi-analytically and introduced directly into the right-hand side of the discrete system.

The remaining domain term $f_i^{\left(n\right)}$ , which contains porosity-reaction coupling and nonlinear residuals, may be handled in several ways. In the present work, a pragmatic incremental treatment is adopted. Because the nonlinear evolution is localized and varies gradually within a single time step, $f_i^{\left(n\right)}$ is evaluated from the previous time level and represented either as a piecewise-constant internal field over a small set of subcells or transformed into equivalent concentrated source terms at selected internal points. This preserves the boundary-dominant character of the method while avoiding full-domain meshing. Such a strategy is especially appropriate for reactive carbonate injection, where strong property evolution is typically confined to regions near the injector and along preferential flow paths.

3.5. Subdomain decomposition and interface integral equations

The reservoir is partitioned into $N_r$ subregions,

For each subregion ${\mathrm{\Omega }}_i$ , a separate boundary integral equation is written over its boundary

where ${\mathrm{\Gamma }}_i^{\mathrm{ext}}$ is the external boundary segment and ${\mathrm{\Gamma }}_{\mathit{int}}^{\mathit{int}}$ denotes interfaces shared with adjacent regions. Across each interface ${\mathrm{\Gamma }}_{\mathit{ij}}=\partial {\mathrm{\Omega }}_i\cap \partial {\mathrm{\Omega }}_j$ , continuity of pressure and flux is imposed:

These relations ensure exact hydraulic transmission across region boundaries. They also play a decisive role in assembling a single global algebraic system from the local subdomain equations. Because the method enforces interface conditions directly at the boundary level, heterogeneous reservoir structure can be captured without volumetric matching grids. This feature is one of the key strengths of BEM for reservoir applications with strong material contrasts.

3.6. Temporal discretization and history integration

To advance the solution in time, the interval $[0,T]$ is divided into $N_t$ steps, and the boundary variables are approximated over each increment using a suitable interpolation rule. For robustness and simplicity, the present formulation adopts a backward time-marching scheme, with boundary pressure and flux assumed constant or linearly varying over each step.

Denote by $p_i^{n+1}\text{and}{q}_{n,i}^{n+1}$ the boundary unknowns at the end of the $(n+1)-\mathrm{th}$ time step. The time convolution integrals involving $G_i\text{and}\partial G_i/\partial n$ are then evaluated over the current step and, when needed, accumulated from previous steps through a history term. The generic form becomes

More explicitly,

3.7. Spatial discretization of the boundary

Each boundary segment Γi\Gamma_iΓi is discretized into NiN_iNi boundary elements. For clarity and computational efficiency, constant or linear interpolation may be adopted. In the constant-element formulation, the boundary pressure and normal flux are assumed uniform over each element $e:$

Substituting these approximations into the boundary integral equation and collocating at the element centers yields

For non-singular element interactions, these integrals are evaluated by standard Gaussian quadrature. For singular or near-singular integrals, analytical evaluation or singularity subtraction is employed to maintain numerical accuracy. Because only boundaries and interfaces are discretized, the resulting system is significantly smaller than the corresponding finite element or finite difference system for the same physical domain.

If higher accuracy is required near wells, corners, or strong permeability gradients, the present constant-element approximation can be replaced by discontinuous linear or quadratic elements without modifying the general structure of the formulation.

3.8. Matrix form of the subdomain equations

Collecting all collocation equations for a subregion Ωi\Omega_iΩi , the discrete system may be written compactly as

The right-hand-side vector ${\mathbf{b}}_i^{(n+1)}$ contains:

Because the matrices are assembled on a subdomain basis, the method naturally accommodates spatial heterogeneity and localized coefficient evolution. This is particularly important for carbonate reservoirs, where dissolution-driven permeability changes may be concentrated in only a limited portion of the domain.

3.9. Global assembly with interface continuity constraints

The global reservoir system is obtained by assembling all subdomain equations and imposing interface continuity conditions. Let

After merging common interface unknowns and enforcing

the assembled system assumes the block form

This fully coupled matrix system is solved at each time step to recover the pressure field on all boundaries. Once the boundary solution is known, interior pressures may be evaluated at arbitrary observation points by substituting the solved boundary data back into the integral representation formula. This post-processing step requires no additional mesh generation and provides a continuous semi-analytical field reconstruction, which is another important advantage of the BEM framework.

3.10. Nonlinear update algorithm for reaction, porosity, and permeability

Because the boundary matrices depend on coefficients that evolve with the reactive state, the solution proceeds incrementally in time. At the end of each time step, the pressure solution is used to update the reaction progress, porosity, and permeability. A staggered algorithm is adopted as follows.

3.10.1. Step 1: solve the linearized BEM system

Given ${\xi }^n,{\varphi }^n,\text{and}{k}^n,$ compute the effective storage and diffusivity in each subregion and solve

3.10.2. Step 2: update the reaction progress variable

Using the pressure field obtained in Step 1, update $\xi$ according to the reaction law introduced in Section 2. For example, with an implicit or semi-implicit Euler rule,

This update preserves the irreversible nature of carbonate dissolution and is straightforward to implement.

3.10.3. Step 3: update porosity and permeability

The porosity is then updated from

or, in the reduced model,

The permeability is updated using

or, in the reduced reactive-flow model,

3.10.4. Step 4: convergence check within the time step

If strong nonlinear coupling is present, the coefficient update may be iterated within the same time step until

This incremental strategy is consistent with the problem statement in the manuscript, where nonlinear permeability evolution is captured through time-stepwise linearization while maintaining the computational efficiency of the boundary element framework.

3.11. Numerical characteristics of the proposed BEM formulation

The present discretization possesses several properties that are particularly advantageous for reactive CO₂ injection modeling.

First, the method remains predominantly boundary-only, thereby avoiding the full volumetric meshing burden of finite difference and finite element approaches. This is especially beneficial for large, laterally extensive reservoirs with complex internal interfaces. Second, the use of subdomain decomposition allows sharp heterogeneity contrasts to be captured rigorously through interface continuity equations. Third, internal wells are represented analytically as singular source terms, eliminating the need for local grid refinement around injection points. Fourth, the nonlinear permeability evolution is incorporated through an incremental update procedure rather than repeated remeshing, which substantially enhances computational efficiency.

At the same time, the method retains the analytical rigor of Green’s-function-based transient diffusion solutions. This combination of dimensional reduction, exact interface treatment, and nonlinear coefficient updating is what enables the proposed framework to model pressure diffusion, reaction-driven permeability enhancement, and evolving injectivity in heterogeneous carbonate reservoirs with high efficiency and accuracy. These capabilities are fully consistent with the objectives and contributions stated in the manuscript.

3.12. Reduced discrete form used in the present implementation

For the numerical examples presented later in the paper, the full formulation is implemented in reduced form by neglecting explicit mechanical equilibrium within the BEM solve and incorporating geomechanical influence through the updated permeability law. Accordingly, the pressure subproblem solved at each step is

The corresponding BEM system at time level $n+1$ is written as

Here, ${\mathbf{b}}_{\text{react}}^{(n+1)}$ contains the equivalent forcing associated with reaction-induced storage change, while the coefficient matrices are assembled using the current subdomain diffusivities. This reduced implementation preserves the essential coupling between pressure diffusion and reaction-driven permeability evolution and forms the computational core of the present study.

4.1. Mathematical validity of the boundary integral representation

The starting point for validity is the observation that the transient boundary integral equation derived in Section 3 is obtained directly from the linearized diffusion equation through Green’s second identity and the adjoint fundamental solution. Consequently, for each time increment and each subregion ${\mathrm{\Omega }}_i,$ the integral representation is mathematically equivalent to the corresponding linear diffusion subproblem, provided that the solution remains sufficiently regular, and the boundary and interface data satisfy the prescribed continuity conditions.

Within a single increment, the local pressure field satisfies

This point is important from a methodological perspective. The proposed scheme is not an empirical boundary approximation of a nonlinear reservoir model; rather, it is a rigorously derived boundary integral reformulation of a sequence of incrementally linearized governing equations. The mathematical validity of the formulation therefore follows from three facts:

Accordingly, if the incremental update converges, the solution obtained by the proposed BEM procedure is a consistent approximation of the original nonlinear reactive-flow problem.

4.2. Consistency of the incremental linearization

The principal approximation introduced in the present formulation arises from the treatment of nonlinearities through time-stepwise linearization. It is therefore necessary to show that the adopted staggered update remains consistent with the governing equations as the time increment $\mathrm{\Delta }t$ decreases.

Let the nonlinear hydraulic operator be written abstractly as

Over a time step, the operator is approximated by

This result establishes the consistency of the proposed nonlinear treatment. It also clarifies the role of time-step selection. Large time steps may still yield stable solutions, but strong reaction fronts or rapid permeability changes require smaller increments to preserve local accuracy of the coefficient update. Thus, the validity of the approach is not merely formal; it is tied to a numerically controlled approximation whose error can be reduced systematically through temporal refinement.

4.3. Recovery of limiting benchmark solutions

A strong verification argument is obtained by demonstrating that the proposed formulation reduces correctly to known benchmark problems under special parameter choices. The following limiting cases are especially relevant.

4.3.1. Pure pressure diffusion in a homogeneous reservoir

If the reaction is suppressed by setting $\lambda =0,$ then

In this case, the proposed BEM becomes identical to a standard transient reservoir boundary element formulation and should recover the well-known radial pressure response in infinite or bounded domains. Agreement in this limit verifies the correctness of the fundamental solution, time convolution treatment, and singular-source representation.

4.3.2. Heterogeneous but nonreactive subdomain diffusion

If $\lambda =0\mathrm{but}{k}_0\left(\mathbf{x}\right)$ varies between subregions, the formulation reduces to transient diffusion in a piecewise heterogeneous reservoir. The exact enforcement of

4.3.3. Spatially uniform reactive evolution

If the reservoir is homogeneous and the pressure field varies slowly enough that the reaction term becomes effectively uniform over a region, then the permeability evolution law

The ability of the proposed formulation to recover these limiting solutions is a fundamental indicator of correctness, because it shows that the method behaves appropriately when the full complexity of the reactive–heterogeneous system is reduced to classical special cases. These limiting cases establish the first level of validation by showing that, in the absence of reactive coefficient evolution, the proposed formulation reduces correctly to classical transient diffusion behavior in both homogeneous and heterogeneous settings.

4.4. Verification against analytical and semi-analytical reference solutions

For practical verification, the first level of quantitative testing should be performed against analytical or semi-analytical solutions in idealized geometries. The most natural benchmark is the transient pressure response to constant-rate injection in an infinite homogeneous reservoir. In the absence of reaction, the classical logarithmic or exponential-integral pressure solution provides a direct reference for evaluating the BEM pressure field at selected observation points and times.

Let $p^{\mathrm{ref}}(\mathbf{x},t)$ denote the reference solution and $p^{\mathrm{BEM}}(\mathbf{x},t)$ the numerical result. The pointwise relative error may be defined as

For weakly reactive cases, semi-analytical reference solutions may be constructed by operator splitting or by solving the one-dimensional reduced problem with a highly refined finite difference grid. In that setting, the proposed BEM solution should reproduce both the pressure profile and the monotonic increase of permeability near the injection region. Agreement under such conditions provides evidence that the incremental nonlinear update remains accurate in regimes where no exact closed-form solution is available.

Having established benchmark consistency in the nonreactive limit, the next step is to verify that the numerical approximation converges systematically under boundary and time-step refinement before assessing agreement with independent reference methods and published reactive data.

4.5. Spatial convergence of the boundary discretization

After benchmark consistency is established, the second level of validation concerns numerical convergence, namely whether the boundary and time discretizations can be systematically refined to approach a stable solution.

The accuracy of the proposed method also depends on the spatial approximation of boundary variables. Because the interior domain is represented analytically through Green’s functions, the principal spatial error arises from discretization of the external boundary, internal interfaces, and any auxiliary internal source representation used for nonlinear correction terms.

Let $N_b$ denote the total number of boundary elements. For a sufficiently smooth boundary solution, the discrete approximation is expected to converge as

A practical mesh-convergence study may be performed by successively refining the boundary and interface discretization and monitoring a scalar response quantity such as:

Let $J_h$ denote one such quantity computed on a mesh of characteristic boundary size $h$ . The relative discretization error between two successive meshes may be estimated as

Convergence is established when ${\eta }_h$ decreases monotonically with refinement and the solution approaches a mesh-independent value. In reactive cases, convergence must be assessed for both hydraulic and constitutive quantities, especially $p,\xi ,\text{and}k$ , because small pressure errors may become amplified through exponential permeability updates.

Accordingly, spatial convergence is treated here as the first numerical proof that the boundary discretization does not merely produce plausible pressure fields, but approaches a mesh-independent solution in a controlled manner.

4.6. Temporal convergence and stability of the time-marching scheme

Because the present formulation employs transient Green’s functions and incremental coefficient updating, temporal accuracy plays a decisive role in overall solution quality. For the adopted backward or semi-implicit time integration, the local truncation error is of first order in time, and therefore one expects

A temporal refinement study is performed by computing the solution with decreasing time steps $\mathrm{\Delta }t,\frac{\mathrm{\Delta }t}{2},\frac{\mathrm{\Delta }t}{4},$ and measuring the change in selected outputs such as well pressure, reaction progress near the injector, or pressure-front location. If the method is consistent, the differences between successive temporal approximations should decrease in accordance with the expected order of accuracy.

Stability is assessed by examining whether the numerical solution remains bounded and physically admissible over long injection periods. In particular, the scheme should satisfy the following qualitative properties:

These constraints are not merely physical requirements; they also serve as indicators that the staggered nonlinear updates have not introduced spurious oscillations or instability. In practice, unconditional linear stability of the implicit hydraulic step is an important advantage of the present approach, while the nonlinear reaction update may require adaptive time-step control when permeability evolves rapidly.

Temporal convergence is equally important in the present reactive problem because time-stepping errors affect not only pressure history but also the subsequent constitutive update of porosity and permeability.

Together, the spatial and temporal refinement requirements define the second level of validation by showing that the proposed approximation is systematically controllable and not dependent on an arbitrary discretization choice.

4.7. Flux conservation and interface accuracy

For heterogeneous subdomain problems, one of the most sensitive indicators of numerical accuracy is the correct transmission of hydraulic flux across interfaces. Because the BEM formulation enforces continuity conditions directly at internal boundaries, the discrete solution should satisfy pressure continuity and local flux balance to within numerical quadrature and algebraic solution tolerance.

For any interface ${\mathrm{\Gamma }}_{\mathit{ij}}$ , the pressure mismatch may be quantified by

A valid and accurate assembly procedure should drive both quantities toward zero under mesh refinement. In addition, global mass conservation over the entire computational domain should be checked by verifying that

Even in the reduced incompressible-density approximation, an analogous storage-plus-flux balance should hold numerically. This global check is particularly valuable in reactive simulations because constitutive changes in porosity may create apparent mass imbalance if implemented inconsistently.

4.8. Accuracy of the singular-source representation

The injection well is modeled as an internal singular source through the fundamental solution. This treatment is one of the key efficiencies of the proposed BEM framework, but it also requires verification to ensure that near-well pressure behavior is captured accurately.

The validity of the singular-source representation may be assessed by comparing the computed pressure field near the injection point with:

In particular, if the pressure field is evaluated at a sequence of points approaching the well, the numerical solution should reproduce the expected radial singular structure without spurious oscillation. Because BEM embeds the singularity analytically through the Green’s function, it is especially well suited to this task. Accurate recovery of near-well pressure gradients is important not only for hydraulic fidelity but also because reaction progress and permeability enhancement are typically strongest in that region.

4.9. Influence of nonlinear permeability evolution on numerical accuracy

The most demanding aspect of the present formulation is the nonlinear feedback between pressure, dissolution, and permeability. Even when the linearized diffusion problem is solved accurately within a time step, errors in pressure may be propagated into the constitutive update through

Because of the exponential dependence, moderate pressure or reaction errors may translate into larger permeability discrepancies over long times. The accuracy of the full coupled solution must therefore be evaluated not only in terms of pressure but also in terms of constitutive quantities.

A useful nonlinear error indicator is the relative coefficient update

If $R_k^{(m+1)}$ decreases monotonically and falls below a prescribed tolerance, the staggered update may be regarded as converged for that time step. In practical terms, this provides a direct diagnostic of whether the chosen time increment is sufficiently small for the reactive front dynamics being simulated.

4.10. Comparison with domain-based numerical methods

An important element of validation is comparison with an established domain discretization technique, such as the Finite Element Method or Finite Difference Method. Such a comparison serves two purposes. First, it provides an external benchmark when analytical solutions are unavailable. Second, it demonstrates that the proposed boundary-only approach reproduces the same physical solution while requiring substantially fewer spatial degrees of freedom for large reservoir domains.

For a representative heterogeneous injection problem, one may compare:

If the BEM and a sufficiently refined domain-based method yield close agreement in pressure and reaction fields, then the proposed formulation may be regarded as accurate from an engineering standpoint. Differences, when present, are most likely to appear near strongly localized nonlinear zones, where the treatment of equivalent domain terms becomes decisive. Such comparisons are therefore especially useful for calibrating the internal source representation adopted for the nonlinear residual.

In the present study, this reference-method comparison is summarized quantitatively in Table 3, which shows that the proposed BEM reproduces the benchmark pressure response with lower maximum error than the compared finite-difference and finite-element implementations while requiring substantially fewer computational resources.

This cross-method agreement provides the third level of validation by demonstrating that the proposed boundary-based formulation reproduces the same physical response as established domain-based methods in a representative storage benchmark.

4.11. Error budget and numerical sensitivity of the proposed BEM formulation

To complement the consistency and benchmark arguments developed in Sections 4.5–4.10, it is useful to identify the dominant numerical error sources in the present implementation and to relate them directly to the observed validation results. In the proposed BEM framework, the total numerical error does not arise from a single approximation, but from the combined effect of boundary interpolation, temporal discretization, staggered nonlinear updating, approximation of equivalent domain terms, and quadrature treatment of regular and singular integrals. The relative importance of these contributions depends on the simulation regime. In nonreactive benchmark cases, the error is governed primarily by boundary resolution and singular-source evaluation, whereas in reactive simulations the dominant contribution progressively shifts toward time integration and nonlinear constitutive updating because permeability depends exponentially on reaction progress and pressure. This distinction is consistent with the validation strategy established in Sections 4.5–4.10 and with the reactive trends observed later in Figures 8–10 and Tables 5–6.

Figure 1. Conceptual model of reactive CO₂ injection in a heterogeneous carbonate reservoir.

Figure 1 presents the physical problem addressed in the study. A CO₂ injection well is located in a heterogeneous carbonate reservoir containing regions of contrasting permeability. Pressure diffuses outward from the well, while CO₂–rock interaction creates a reactive zone in which dissolution modifies pore structure and increases permeability. The physical meaning of the figure is that reservoir behavior is controlled by a coupled feedback loop: injection raises pressure, pressure promotes reaction, reaction alters permeability, and the evolving permeability field changes subsequent flow and injectivity.

Figure 2. Computational workflow of the proposed reactive boundary element formulation.

Figure 2 summarizes the numerical solution strategy used in the paper. The workflow starts from reservoir inputs and governing equations, proceeds through transient BEM solution of the pressure field, and then updates reaction progress, porosity, and permeability at each time step. The physical meaning is that the reservoir does not evolve under fixed properties; instead, hydraulic and reactive processes are solved in a stepwise coupled manner so that pressure and material-property changes influence each other throughout the simulation.

Figure 3. Temporal evolution of carbonate dissolution reaction.

Figure 3 shows the growth of the reaction progress variable with time and the corresponding decrease in reaction rate. Initially, the reaction proceeds rapidly because the injected CO₂ creates strong disequilibrium, but the rate slows as the system evolves. The physical meaning is that carbonate dissolution is strongest during the early stage of injection and gradually approaches saturation, which is consistent with irreversible but self-limiting reactive behavior.

Figure 4. Reaction-induced permeability enhancement.

Figure 4 shows how normalized permeability increases as reaction progress advances for different values of the permeability-growth coefficient. The physical meaning is that dissolution enlarges the pore network and improves hydraulic conductivity, but the strength of this effect depends on how sensitive permeability is to chemical alteration. This relationship is central to the model because it links geochemical change directly to injectivity and pressure evolution.

Figure 5. Schematic heterogeneous reservoir configurations used in the nonreactive benchmark cases.

Figure 5 shows the two heterogeneous reservoir layouts used to isolate the effect of permeability architecture under nonreactive conditions. One case contains a high-permeability channel intersecting the injector, while the other contains a low-permeability barrier surrounding the near-well region. The physical meaning is that heterogeneity alone can either facilitate or restrict pressure diffusion, even when no dissolution or permeability evolution is present.

Figure 6. Pressure redistribution in heterogeneous nonreactive reservoirs.

Figure 6 presents pressure contours for the two heterogeneous nonreactive cases. In the channel case, pressure propagates preferentially along the conductive pathway, whereas in the barrier case pressure remains concentrated near the well. The physical meaning is that the initial permeability structure controls the spatial pattern of pressure migration and therefore strongly influences early-time injectivity and reservoir pressurization.

Figure 7. Wellbore pressure histories for nonreactive homogeneous and heterogeneous cases.

Figure 7 compares the temporal evolution of wellbore pressure buildup for the homogeneous reference case, the high-permeability channel case, and the low-permeability barrier case. The channelized case yields the smallest pressure buildup, while the barrier case yields the largest. These results show that static heterogeneity directly controls pressure dissipation efficiency.

Figure 8. Reactive evolution in a homogeneous carbonate reservoir.

Figure 8 shows the spatial distributions of pressure, reaction progress, and normalized permeability in the reactive homogeneous case. A localized reactive halo forms around the injector, where pressure is highest, and dissolution is most intense. The physical meaning is that even in a geologically uniform reservoir, reactive CO₂ injection creates a nonuniform altered zone near the well, and this localized change governs the main hydraulic feedback.

Figure 9. Pressure response in homogeneous nonreactive and reactive cases.

Figure 9 compares wellbore pressure buildup in the homogeneous reservoir with and without reactivity. At early time the two cases are similar, but later the reactive case shows lower pressure buildup. The physical meaning is that dissolution-driven permeability enhancement reduces hydraulic resistance near the injector and therefore moderates long-term pressurization.

Figure 10. Injectivity evolution in homogeneous cases.

Figure 10 shows the apparent injectivity index for the nonreactive and reactive homogeneous cases. Injectivity increases more strongly in the reactive case as permeability evolves near the well. The physical meaning is that geochemical alteration can improve injection performance over time by making the near-well region more conductive.

4.11.1. Boundary interpolation error

Boundary interpolation error arises from approximating pressure and normal flux on the external boundary and internal interfaces by constant or low-order element fields. In the present formulation, this source of error is directly linked to the spatial convergence discussion of Section 4.5, where the principal discretization error is attributed to the boundary and interface representation rather than to volumetric meshing. The numerical evidence reported in Table 3 shows that, with only 320 boundary elements, the proposed BEM attains a maximum pressure error of 0.8%, compared with 1.6% for the finite-difference model using a 200 × 200 grid and 1.2% for the finite-element model using 160,000 elements. This comparison indicates that the boundary interpolation adopted in the present study is already sufficiently accurate for the benchmark configuration and that the boundary-only discretization does not compromise pressure accuracy.

In heterogeneous cases, the adequacy of the boundary interpolation is further supported by the absence of spurious field discontinuities across internal interfaces. The manuscript already notes that pressure continuity and flux continuity are preserved across strongly contrasting subdomains without visible jumps in the reconstructed solution, which indicates that interpolation error at interfaces remains controlled even when the permeability architecture is nonuniform.

4.11.2. Temporal discretization error

Temporal discretization error arises from approximating the transient convolution integrals and from advancing the reaction and property fields over finite time increments. As discussed in Section 4.6, this contribution is expected to be first order in time for the adopted backward or semi-implicit scheme, and its practical influence becomes more pronounced when the constitutive variables evolve rapidly. In the present simulations, the importance of temporal resolution is reflected by the cumulative separation between the nonreactive and reactive responses in Table 6: the pressure reduction associated with reactive permeability evolution increases from 2.1% at 10 days to 18.1% at 120 days, while the injectivity gain increases from 2.9% to 21.4% over the same interval. These trends confirm that even modest time-integration errors in the early pressure field may accumulate through the coupled update and become more significant at later times.

Accordingly, temporal discretization is not merely a secondary implementation detail in the present reactive problem; it directly affects the predicted onset and growth rate of the altered near-well zone. This interpretation is consistent with the monotonic time evolution reported for reaction progress and permeability in Table 5 and with the stability-oriented discussion given in Section 4.6.

4.11.3. Nonlinear update error

Nonlinear update error is introduced by freezing storage and permeability within each time increment and by terminating the internal staggered iterations at a finite convergence tolerance. Among all error sources, this is the one most closely tied to the distinctive physics of the present model, because the hydraulic coefficients are not fixed but evolve in response to pressure-activated dissolution. Section 4.9 already identifies this issue as the most demanding aspect of the formulation, emphasizing that moderate pressure errors may be amplified through the exponential permeability law.

The magnitude of this sensitivity is visible in the near-well constitutive evolution reported in Table 5. Over the simulated interval, the reaction progress variable increases from 0.12 at 10 days to 0.73 at 120 days, while the normalized permeability increases from 1.18 to 2.40. Because such a permeability change is produced through the staggered constitutive update, the accuracy of the nonlinear iteration directly controls the accuracy of the long-time pressure and injectivity predictions. This is also evident in Table 6, where the reactive pressure response departs progressively from the nonreactive baseline as the nonlinear feedback strengthens.

For this reason, nonlinear update error should be regarded as the principal accuracy limiter in strongly reactive simulations, even when the underlying linear BEM solve is itself highly accurate. This conclusion follows directly from the constitutive trends already documented in Sections 4.9 and 5.6.

4.11.4. Domain-term approximation error

Domain-term approximation error arises because the present BEM formulation represents nonlinear residuals and porosity-related internal effects through equivalent internal source terms rather than through exact full-domain integration. This approximation is methodologically important because it enables the formulation to remain predominantly boundary-only, but it also introduces a potential source of model discrepancy in regions where the reactive front becomes highly localized. Section 4.10 already points out that the main differences between the BEM and domain-based reference methods are expected to appear precisely in such strongly localized nonlinear zones.

The best available quantitative evidence that this approximation remains acceptable is provided by the validation comparisons already incorporated into the manuscript. First, Table 3 shows that the global pressure prediction remains highly accurate despite the reduced domain treatment, with the BEM giving a 0.8% maximum pressure error while remaining markedly cheaper computationally than the finite-difference and finite-element alternatives. Second, Table 1 reports small deviations between the BEM permeability predictions and published carbonate-dissolution reference studies, with reported relative differences of 2.6%, 2.1%, and 1.6% for the three comparison cases. Taken together, these results suggest that the adopted equivalent-source approximation is sufficiently accurate for the class of reactive CO₂ injection problems considered here, although it is likely to become more critical in simulations involving sharper fronts or stronger local channelization.

4.11.5. Quadrature and singular-integration error

Quadrature and singular-integration error arises during evaluation of regular, near-singular, and singular boundary integrals, especially in the vicinity of the injection well and at closely spaced interfaces. This source of error is particularly relevant because the present method represents the well analytically as an internal singular source, so accurate integration is essential to recover the correct near-well pressure gradient. Section 4.8 explicitly identifies the singular-source treatment as one of the key efficiencies of the formulation and emphasizes that near-source behavior must be captured without spurious oscillation.

The quality of the present singular treatment is supported indirectly by two results already reported in the manuscript. In the nonreactive benchmark case, the pressure field remains smooth and physically admissible while reproducing the expected radial diffusion behavior near the injector. In addition, Table 3 shows that the proposed BEM achieves the lowest maximum pressure error among the three compared methods, which would not be possible if singular or near-singular integration errors were dominant in the implementation. The heterogeneous case studies provide further support, since interface transmission is reported to occur without spurious jumps in the reconstructed field.

4.11.6. Relative importance of the error sources

For the present implementation, the error sources do not contribute equally. In the baseline nonreactive benchmarks, the dominant contributions are boundary interpolation and singular-integration accuracy, whereas the favorable results in Table 3 indicate that both are already well controlled at the adopted discretization level. In contrast, once reactive coupling is activated, the dominant error contribution shifts toward the constitutive side of the problem. The progressive permeability increase reported in Table 5 and the growing pressure and injectivity separation reported in Table 6 show that nonlinear update error is the most critical source in the present formulation, followed by temporal discretization error, then domain-term approximation error, while boundary interpolation and quadrature/singular-integration errors are secondary for the benchmark settings considered here. This ranking is consistent with the formulation itself: the BEM operator is accurate and efficient for the linearized diffusion step, but the strongest numerical sensitivity enters through time-dependent reactive property evolution.

4.12. Practical validation criteria for the present study

In the context of the present paper, the proposed BEM formulation is considered valid and accurate when the following criteria are simultaneously satisfied:

These criteria collectively address both the mathematical and engineering dimensions of validity. They ensure that the method is not only formally derived from the governing equations, but also capable of delivering quantitatively reliable predictions for injection-induced pressure buildup and permeability evolution in heterogeneous carbonate formations.

The fourth level of validation addresses the reactive component of the model, since storage-relevant predictive value depends not only on hydraulic accuracy but also on whether reaction-induced permeability evolution remains physically credible.

For this reason, Table 1 compares the predicted permeability enhancement with representative published carbonate-dissolution studies and serves as the literature-based reactive validation step of the present framework.

The benchmark comparisons, interface-consistency observations, constitutive trends, and cross-method agreement presented in Tables 1 and 3 and in Sections 4.5–4.11 indicate that these criteria are satisfied to a level sufficient for the present reduced reactive-flow implementation.

Agreement with these published trends provides the fourth level of validation by showing that the reactive constitutive response of the model is consistent with established carbonate-dissolution behavior relevant to CO₂ storage.

4.13. Concluding assessment of validity and accuracy

The preceding analysis establishes the validity and numerical reliability of the proposed formulation through four complementary levels of evidence: benchmark recovery, convergence requirements, agreement with independent reference methods, and literature-based validation of reactive permeability evolution. At the continuous level, the method follows rigorously from the time-dependent Green’s-function representation of the incrementally linearized diffusion problem. At the discrete level, accuracy is controlled by boundary refinement, time-step selection, interface enforcement, and nonlinear update tolerance. The formulation correctly recovers classical nonreactive limits, accommodates heterogeneous subdomains through exact interface conditions, and incorporates evolving permeability through a consistent staggered update strategy.

Accordingly, the proposed method provides a reliable computational framework for the simulation of reactive CO₂ injection in carbonate reservoirs, provided that the benchmark checks and convergence criteria outlined above are satisfied in implementation. This establishes a firm foundation for the numerical experiments presented in the following section, where the formulation is applied to representative heterogeneous storage scenarios in order to examine pressure diffusion, reaction-induced permeability enhancement, and their implications for storage performance.

Taken together, the results of this section demonstrate that the proposed formulation satisfies the main requirements of a reliable storage-performance model for reactive CO₂ injection. It recovers the correct limiting benchmark behavior, admits systematic spatial and temporal refinement, agrees with independent reference methods, and reproduces reactive permeability trends consistent with published carbonate-dissolution studies. The formulation can therefore be considered sufficiently validated for the class of heterogeneous reactive storage problems addressed in this work.

5.1. Computational model and reservoir configuration

All numerical examples are based on a two-dimensional horizontal carbonate reservoir containing a centrally located injection well represented analytically as an internal singular source. This representation preserves the principal advantage of the BEM, namely boundary-only discretization, while retaining correct near-well pressure behavior. The reservoir is assumed to be fully saturated and laterally extensive enough that the external boundary does not dominate the early-time local response. The initial state is defined by a uniform reference pressure, while the reservoir may be homogeneous or partitioned into multiple subregions with distinct initial permeability values depending on the scenario considered. In the reduced implementation used here, the governing model consists of the transient hydraulic diffusion equation coupled to the pressure-activated reaction law and the constitutive update laws for porosity and permeability. The parameter set adopted for the simulations is consistent with the physically realistic assumptions summarized before Figure 1, including the moderate kinetic rate, the finite simulation horizon, and the exponential permeability-growth law appropriate for carbonate dissolution.

For most examples, a far-field pressure condition is imposed on the outer boundary in order to approximate hydraulic communication with an undisturbed surrounding formation, although alternative no-flow segments are introduced in selected cases to assess confinement effects. The monitoring strategy adopted in the present section is directly connected to the earlier response metrics illustrated in Figures 1–5, where pressure evolution at the injection point, delayed pressure response at a distant observation point, reaction progress, permeability enhancement, and radial pressure distribution were first established as the key physical indicators of model behavior. These earlier figures therefore serve not merely as preliminary illustrations, but as reference diagnostics for interpreting the more complex scenarios considered below.

5.2. Discretization strategy and numerical implementation

The external boundary of the reservoir and all internal material interfaces are discretized using constant boundary elements, and the transient response is obtained through the implicit time-marching procedure derived in Section 3. At each time step, the BEM matrices are assembled using the current subdomain diffusivities, the boundary pressure and flux are solved, the interior pressure field is reconstructed at selected monitoring points, the reaction progress variable is updated, and the porosity and permeability are then revised through the constitutive laws. When nonlinear coupling becomes pronounced, a small number of internal corrective iterations is performed within the same time increment until the coefficient update satisfies the prescribed convergence criterion. This implementation follows the staggered update framework described in the draft and preserves the analytical structure of the boundary-integral problem while accommodating time-dependent hydraulic properties.

Unless otherwise stated, the baseline simulations reported below use the boundary discretization and time-stepping configuration associated with the benchmark setting summarized in Table 3.

The robustness of this discretization strategy is already foreshadowed by the smooth pressure traces shown in Figures 1 and 2 and by the monotonic constitutive evolution reported in Figures 3 and 4. Moreover, the favorable comparison with finite difference and finite element methods in Table 3 indicates that the proposed BEM achieves the required level of accuracy with substantially fewer unknowns and a markedly lower computational cost. The same numerical framework is therefore adopted consistently in all subsequent case studies and sensitivity analyses.

5.3. Definition of response indicators

To facilitate comparison between scenarios, four response measures are used throughout this section. The first is the wellbore pressure buildup, defined as the pressure increment at or near the injector relative to the initial state; this quantity provides the most direct indicator of injectivity limitation and pressurization risk. The second is the reaction-zone extent, defined as the spatial region in which the reaction progress exceeds a prescribed threshold; its corresponding area is used as a measure of the growth of the chemically altered zone. The third is the mean permeability enhancement, evaluated either globally or within a near-well subregion, and used to quantify the hydraulic amplification produced by dissolution. The fourth is the apparent injectivity index, defined as the ratio of injection rate to wellbore pressure buildup, which measures the degree to which local hydraulic performance improves during sustained injection. These indicators are consistent with the validation trends established in Tables 1 and 3 and with the performance comparisons later reported in Tables 4–8.

5.4. Example 1: Baseline nonreactive homogeneous reservoir

The first example considers a homogeneous, nonreactive reservoir and serves as the reference case for the entire numerical study. The initial permeability and porosity are spatially uniform, and the reactive contribution is suppressed so that the problem reduces to the classical transient diffusion equation with constant diffusivity. Under these conditions, the numerical solution provides a direct test of the transient BEM implementation and establishes the baseline against which the effects of heterogeneity and reaction can be measured.

The resulting pressure response exhibits the expected radially symmetric diffusion behavior centered on the injection well. The near-well pressure increase follows the monotonic logarithmic-type trend shown in Figure 1, confirming correct representation of point-source injection and transient pressure buildup. At a more distant observation point, the delayed arrival of the diffusion front is consistent with the response displayed in Figure 2. The spatial pressure field remains smooth and physically admissible, and the radial pressure decay at a representative later time is consistent with the profile shown in Figure 5. Together, Figures 1, 2, and 5 confirm that the nonreactive solution captures the correct early-time near-well buildup, delayed far-field propagation, and spatial decay of pressure. The computational efficiency of the BEM in this baseline case is summarized in Table 3, which shows that the method achieves accuracy comparable to or better than finite difference and finite element reference models at substantially lower computational cost.

5.5. Example 2: Nonreactive heterogeneous reservoir

The second example introduces geological heterogeneity while keeping the system nonreactive so that the role of the initial permeability architecture can be isolated. The reservoir is partitioned into multiple subregions with distinct permeability values, but the reaction variable remains inactive, and the permeability field does not evolve in time. Two representative configurations are considered: a conductive channel intersecting the injector and a low-permeability barrier or annular damaged zone surrounding the near-well region. These subdomain arrangements correspond to the heterogeneous layouts summarized in Table 4.

The heterogeneous pressure fields show a strong departure from the radial symmetry observed in the baseline case. In the channelized configuration, pressure propagates preferentially along the more conductive path, producing elongated contours and a reduced wellbore pressure buildup; this type of behavior is represented in Figure 6a. In contrast, when the injector is partially enclosed by a low-permeability barrier, the pressure field becomes more localized, gradients steepen near the interfaces, and the wellbore pressure rises more rapidly, as shown schematically in Figure 6b. The corresponding pressure histories are compared in Figure 7, which demonstrates that the channelized case yields the smallest pressurization, whereas the barrier-controlled case produces the largest buildup. The associated quantitative differences in pressure response and subdomain behavior are summarized in Table 4.

This example also verifies the quality of the subdomain BEM assembly. Pressure continuity and flux continuity are preserved across the internal interfaces without spurious jumps in the reconstructed field, confirming that the interface treatment remains accurate under strong permeability contrast. More importantly, the results establish that static heterogeneity alone can substantially alter pressure migration and injection-induced pressurization, even before reactive alteration becomes active.

5.6. Example 3: Reactive homogeneous reservoir

The third example activates the reaction model in an initially homogeneous reservoir in order to isolate the effect of dissolution-driven property evolution. The initial permeability remains spatially uniform, but the reaction progress variable is now allowed to evolve according to the pressure-activated kinetic law. At early times, the solution remains close to the nonreactive baseline because the reaction progress is still small, and the hydraulic coefficients have not yet changed appreciably. As injection continues, however, the elevated near-well pressure activates dissolution, causing the reaction progress variable to increase, followed by local porosity enlargement and permeability enhancement.

The temporal evolution of the reaction variable is represented by Figure 3, which shows rapid early-time growth followed by asymptotic saturation, while Figure 4 shows the corresponding monotonic increase in normalized permeability. The coupled spatial distributions of pressure, reaction progress, and normalized permeability in the reactive homogeneous case are illustrated in Figure 8, and the near-well average values of these quantities are summarized in Table 5. A key feature of this case is the emergence of a reactive halo around the injector, within which the reaction proceeds most rapidly and the permeability grows most strongly. As a consequence, the pressure-history curve departs progressively from the nonreactive baseline, as demonstrated in Figure 9. The associated increase in apparent injectivity is shown in Figure 10, while the quantitative changes in wellbore pressure buildup and injectivity are reported in Table 6.

These results demonstrate that dissolution does not merely add a chemical correction to a hydraulic problem; rather, it modifies the near-well hydraulic architecture in a dynamically coupled manner. Pressure buildup accelerates reaction, reaction increases permeability, and the enhanced permeability in turn moderates subsequent pressure growth. This pressure–reaction–permeability feedback is one of the central mechanisms captured by the present formulation.

5.7. Example 4: Reactive heterogeneous carbonate reservoir

The fourth example combines geological heterogeneity and reactive dissolution and therefore represents the most realistic scenario considered in the present study. The initial permeability is piecewise heterogeneous, and the spatially variable reaction field modifies that architecture progressively as injection proceeds. The heterogeneous reservoir configuration and the corresponding material parameters are summarized in Table 7, while the initial geometric layout is represented in Figure 11.

Figure 11. Heterogeneous carbonate reservoir used in the reactive case.

Figure 11 presents the representative heterogeneous reservoir geometry used in the reactive simulations, including conductive and restrictive subregions and the reactive near-well zone. The physical meaning is that realistic carbonate storage formations are spatially nonuniform, and this initial structure determines where pressure accumulates and where reactive alteration is most likely to develop.

At early times, the pressure response is governed mainly by the initial permeability field: pressure diffuses rapidly through conductive subregions and remains concentrated near restrictive zones. This asymmetry is evident in the early-time pressure contours of Figure 12a. As injection continues, however, the reaction front develops nonuniformly, producing spatially variable permeability enhancement, as illustrated in Figures 12b and 12c. The wellbore pressure history for the reactive heterogeneous case, compared with the nonreactive heterogeneous and reactive homogeneous responses, is presented in Figure 13. In some subregions, reaction-induced permeability growth partially offsets the effect of an initially restrictive permeability distribution and reduces the long-time pressure buildup. In other regions, pre-existing conductive pathways become even more dominant as dissolution further amplifies transmissivity. The evolution of the apparent injectivity index and altered-zone extent is shown in Figure 14, while the principal performance measures are summarized in Table 8.

Figure 12. Reactive evolution in a heterogeneous carbonate reservoir.

Figure 12 shows the spatial distributions of pressure, reaction progress, and normalized permeability in the reactive heterogeneous case. The altered zone is asymmetric because the initial facies architecture redistributes pressure and therefore controls where dissolution is strongest. The physical meaning is that heterogeneity and reaction do not act independently; instead, the initial hydraulic structure guides the later pattern of geochemical modification.

Figure 13. Pressure histories for representative storage scenarios.

Figure 13 compares wellbore pressure buildup for the nonreactive heterogeneous, reactive homogeneous, and reactive heterogeneous cases. The differences between the curves show how both heterogeneity and reactive permeability evolution influence long-term pressurization. The physical meaning is that storage performance depends on the combined effect of geological structure and dynamic property change, not on either one alone.

Figure 14. Performance indicators for the reactive heterogeneous reservoir.

Figure 14 shows the evolution of the apparent injectivity index and the altered-zone extent in the reactive heterogeneous case. As the reactive zone expands, injectivity improves because permeability increases in hydraulically important regions near the injector. The physical meaning is that the growth of the chemically altered zone is directly tied to operational storage performance, especially pressure management and injection efficiency.

This example makes clear that reactive flow in heterogeneous carbonate reservoirs is strongly path dependent. The initial facies structure controls the early pressure field, the pressure field determines where reaction proceeds most actively, and the reaction then redistributes permeability, thereby modifying subsequent flow. The proposed BEM formulation resolves this sequence naturally through subdomain decomposition and time-stepwise coefficient updating, without requiring full-domain remeshing.

5.8. Sensitivity study I: effect of reaction-rate coefficient

The first parametric study examines the influence of the kinetic reaction-rate coefficient while keeping all other parameters fixed. Larger values of the kinetic coefficient lead to faster growth of the reaction progress variable at a given pressure level and therefore accelerate the onset of permeability enhancement. The resulting trends are consistent with the temporal behavior of the reaction variable shown in Figure 3 and with the resulting permeability amplification shown in Figure 4. In practical terms, increasing the reaction rate causes the system to depart earlier from the nonreactive pressure trajectory, reduces the slope of the pressure-history curve sooner, and expands the altered zone more rapidly. These trends are reflected most clearly in the comparative pressure histories of Figure 9 and in the injectivity trends of Figure 10, while the associated numerical values are consistent with Tables 5 and 6.

However, the simulations also indicate that the influence of the reaction-rate coefficient is not indefinite. Because the adopted reaction law includes a saturation factor, the dissolution process gradually slows as the reaction progress approaches its upper limit. The kinetic coefficient therefore primarily controls the timescale of hydraulic alteration rather than changing the fundamental qualitative nature of the long-time response.

5.9. Sensitivity study II: effect of permeability-growth coefficient

The second sensitivity study focuses on the permeability-growth coefficient, which determines how strongly a given level of reaction progress is translated into hydraulic-property enhancement. When this coefficient is small, even appreciable reaction progress produces only mild permeability increase, so the pressure response remains close to that of the original reservoir architecture. When the coefficient is moderate or large, by contrast, relatively small changes in reaction progress produce substantial near-well conductivity enhancement, leading to earlier deviation from the nonreactive pressure trajectory and a more pronounced increase in injectivity. These effects are directly consistent with the constitutive trends in Figure 4, the reactive field distributions in Figure 8, the pressure-history comparisons in Figure 9, and the injectivity evolution shown in Figure 10. The quantitative consequences of varying this coefficient are likewise reflected in Tables 5 and 6.

From a reservoir-engineering viewpoint, this parameter measures how sensitively the pore network responds to carbonate dissolution. Reservoirs characterized by larger permeability-growth coefficients are therefore more likely to display self-enhancing injectivity during sustained CO₂ injection.

5.10. Sensitivity study III: effect of heterogeneity intensity

The third sensitivity study examines the role of the contrast between high- and low-permeability subregions. As the heterogeneity ratio increases, the reservoir exhibits stronger asymmetry in both the pressure field and the reaction field. In the nonreactive setting, this produces more pronounced channeling or shielding effects, as already demonstrated by Figures 6 and 7 and summarized in Tables 3 and 4. In the reactive setting, the same increase in heterogeneity also changes where dissolution is likely to concentrate, since reaction tends to intensify in zones characterized by either stronger pressure buildup or longer hydraulic residence time. This effect is evident in the heterogeneous reactive layouts of Figures 11–14 and in the corresponding performance indicators of Tables 7 and 8.

The results show that strong heterogeneity may either preserve compartmentalization or progressively erode it through localized permeability enhancement, depending on the arrangement of the subdomains and the magnitude of the reactive driving force. This confirms the necessity of retaining explicit subdomain representation in the numerical model when long-term storage behavior is to be forecast in carbonate systems.

5.11. Comparative discussion of pressure evolution

A key conclusion emerging from all scenarios is the distinction between early-time and long-time pressure behavior. At early times, all cases are governed predominantly by the initial permeability field because reaction has not yet significantly altered the pore structure. This is why the early segments of the reactive pressure curves resemble the corresponding nonreactive responses, as seen by comparison of Figures 1, 7, 9, and 13. At intermediate and later times, however, reaction-induced permeability evolution becomes progressively more important and shifts the system into a modified diffusion regime, one that is evident in the constitutive evolution of Figures 3 and 4, the reactive field distributions of Figures 8 and 12, and the injectivity gains shown in Figures 10 and 14. The corresponding numerical trends are reported in Tables 4, 6, and 8.

Thus, the proposed BEM formulation does not merely capture a sequence of static pressure fields; it resolves a genuine transition from an initially imposed hydraulic architecture to a dynamically altered one, with direct consequences for injectivity and pressure management.

5.12. Implications for injectivity and storage performance

The numerical examples indicate that dissolution-driven permeability enhancement can materially reduce injection-induced pressure buildup, especially in the near-well region. From an operational standpoint, this is important because even moderate gains in local conductivity may permit sustained injection at lower pressure risk. This beneficial effect is most clearly reflected in the difference between the nonreactive and reactive pressure histories in Figures 9 and 13 and in the corresponding increase in apparent injectivity shown in Figures 10 and 14. The same trend is supported quantitatively by Tables 6 and 8.

At the same time, the simulations show that permeability enhancement is strongly localized and highly sensitive to reservoir structure. A favorable outcome is most likely when dissolution improves local connectivity without creating uncontrolled preferential pathways. For this reason, Figures 6–14 and Tables 3–8 collectively suggest that pressure management in carbonate CO₂ storage cannot be assessed reliably from static permeability fields alone; rather, it requires simultaneous consideration of both the initial geological architecture and its reactive evolution during injection.

5.13. Computational performance and suitability of the BEM approach

In addition to the physical insights, the numerical experiments confirm the practical advantages of the proposed BEM formulation. Because only the external boundary and internal interfaces are discretized, the total number of unknowns remains comparatively small even in reservoirs with multiple heterogeneous zones. The analytical treatment of the injection well and the post-processed reconstruction of interior fields further reduce computational overhead. This is already evident from Table 2, where the BEM outperforms finite difference and finite element reference models in CPU time and memory demand while preserving high accuracy. The same computational economy makes the method well suited to the multiple scenario comparisons and parametric analyses summarized later in Tables 3–8.

The results therefore confirm that the method provides a favorable compromise between physical realism and numerical efficiency. This is particularly valuable for problems in which the dominant complexities arise from boundary conditions, internal source behavior, and heterogeneous interfaces rather than from the need for a fully volumetric discretization of the entire reservoir.

5.14. Concluding remarks on the numerical results

The numerical results demonstrate that the proposed Boundary Element formulation successfully reproduces classical transient diffusion in homogeneous nonreactive reservoirs, accurately captures the redistribution of pressure caused by geological heterogeneity and resolves the localized porosity and permeability enhancement induced by carbonate dissolution. Specifically, Figures 1, 2, and 5 establish the correct baseline diffusion response; Figures 3 and 4 confirm the physically consistent temporal evolution of reaction progress and permeability; Figures 6 and 7 demonstrate the hydraulic impact of static heterogeneity; Figures 8–10 show how reactive alteration modifies the homogeneous near-well flow regime and improves injectivity; and Figures 11–14 illustrate the path-dependent interplay between heterogeneity and reactivity in the most realistic reservoir configuration. Table 2 documents the computational efficiency of the BEM, Tables 3 and 4 summarize the heterogeneous nonreactive scenarios, Tables 5 and 6 quantify the homogeneous reactive response, and Tables 7 and 8 synthesize the behavior of the reactive heterogeneous reservoir.

Taken together, these findings confirm that long-term CO₂ injection behavior in carbonate formations cannot be described adequately by static permeability alone. Even when the early-time response is governed primarily by the initial facies architecture, sustained reactive injection may substantially reorganize the effective hydraulic structure of the reservoir. The resulting system is nonlinear, path-dependent, and spatially localized, with direct implications for injectivity forecasting, pressure control, and storage-system design.

5.15. Transition to the overall conclusions of the study

The case studies and parametric analyses presented in this section provide strong support for three broader conclusions. First, reaction-induced permeability evolution is not a secondary correction to the hydraulic problem but may become a controlling mechanism in long-duration injection. Second, explicit representation of reservoir heterogeneity is essential because it influences both early pressure redistribution and the later spatial pattern of chemical alteration. Third, the Boundary Element framework provides a computationally efficient and physically rigorous alternative to conventional domain-based methods for this class of reactive storage problems. These observations lead directly to the concluding section, where the principal methodological contributions, physical implications, and future research directions are summarized.

6.1. Principal contributions of the study

This study has developed a new Boundary Element formulation for transient CO₂ injection in heterogeneous carbonate reservoirs with explicit consideration of reaction-induced porosity and permeability evolution. The proposed model was motivated by the need for computationally efficient yet physically meaningful simulation tools capable of capturing the coupled hydraulic and geochemical processes that govern pressure buildup, injectivity, and storage performance in reactive subsurface systems. By embedding the pressure-diffusion problem in a time-dependent boundary-integral framework and coupling it to constitutive laws for reaction progress, porosity growth, and permeability enhancement, the present work extends classical BEM methodology to a nonlinear class of reservoir problems for which boundary-only formulations remain relatively uncommon.

A central contribution of the study lies in the incremental linearization strategy, which preserves the analytical structure of the BEM while allowing hydraulic properties to evolve in time. The method rigorously enforces interface continuity across heterogeneous subregions, treats wells analytically as internal singular sources, and reconstructs transient interior fields without volumetric remeshing. The physical relevance of these methodological advances is demonstrated throughout Section 5 by the consistent interpretation of Figures 1–14 and the quantitative comparisons reported in Tables 2–8. In particular, Table 2 shows that these advances are achieved without sacrificing computational efficiency.

6.2. Main physical findings

The numerical results show that reaction-driven permeability evolution substantially alters pressure development during CO₂ injection. In homogeneous reactive systems, the near-well dissolution halo shown in Figures 3, 4, and 8 produces moderated long-time pressure buildup and increased injectivity, as confirmed by Figures 9 and 10 and by Tables 5 and 6. In heterogeneous reservoirs, the effect is more complex because the initial facies architecture governs early pressure redistribution, while the subsequent reaction redistributes permeability in a spatially nonuniform manner, as illustrated by Figures 11–14 and summarized in Tables 7 and 8. These results demonstrate that the hydraulic response of carbonate storage reservoirs is intrinsically path-dependent and cannot be described adequately using static permeability fields alone.

Equally important, the comparison between homogeneous and heterogeneous nonreactive cases, represented by Figures 5–7 and Table 4, shows that geological heterogeneity exerts a strong control on pressure propagation even before reactive alteration becomes active. Accordingly, reliable interpretation of field-scale injection behavior requires simultaneous treatment of both the initial geological structure and its reactive modification during operation.

6.3. Significance for CO₂ storage assessment

From the standpoint of storage engineering, the findings suggest that carbonate dissolution may, under favorable conditions, improve injectivity by reducing near-well pressure accumulation. This effect is particularly relevant in pressure-limited reservoirs, where even moderate permeability enhancement may lower operational risk and improve sustained injection performance. The improvement in apparent injectivity quantified in Figures 10 and 14 and Tables 6 and 8 therefore has direct practical significance. At the same time, the results also indicate that dissolution-driven enhancement is not uniformly beneficial: in strongly heterogeneous settings, localized permeability growth may reinforce preferential flow paths and modify the effective hydraulic connectivity of the reservoir. For this reason, Figures 6–14 and Tables 3–8 collectively support the conclusion that pressure management and storage forecasting in carbonate systems must be based on coupled reactive-flow simulation rather than on purely hydraulic models with fixed material properties.

6.4. Methodological strengths of the proposed BEM framework

The study further demonstrates that the Boundary Element framework is especially well suited to reactive CO₂-injection problems in which the dominant numerical challenges arise from large reservoir extent, internal source behavior, and subdomain interfaces. The smooth benchmark responses in Figures 1 and 2, the physically consistent constitutive trends in Figures 3 and 4, the accurate treatment of permeability contrasts in Figures 6 and 7, and the efficient handling of evolving heterogeneous structure in Figures 11–14 all indicate that the formulation preserves both numerical stability and physical realism. The performance comparison in Table 3 confirms that this is achieved with substantially lower computational burden than conventional domain-based discretizations. Thus, the present formulation provides a strong basis for large-scale sensitivity analysis, scenario comparison, and long-term forecasting.

6.5. Limitations of the present study