The fluid mass balance in a deformable porous medium may be written in the form ∂ ∂ t ( ϕ ρ f ) + ∇ · ( ρ f v ) = q m , (1) where ρ f is the fluid density, v is the Darcy velocity, and q m denotes a volumetric source term associated with injection or extraction. For slightly compressible flow, the density variation may be linearized as ρ f = ρ f 0 [ 1 + c f ( p − p 0 ) ] , (2) where ρ f 0 is a reference fluid density, c f is the fluid compressibility, and p 0 is a reference pressure. Neglecting higher-order terms and dividing by ρ f 0 , the mass balance reduces to ϕ c f ∂ p ∂ t + ∂ ϕ ∂ t + ∇ · v = q , (3) in which q = q m / ρ f 0 is the normalized source term.

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 S ∂ p ∂ t + α ϕ ∂ ϕ ∂ t + ∇ · v = q , (4) where S incorporates the combined effects of fluid compressibility, pore compressibility, and matrix compressibility, and α ϕ is a coupling coefficient associated with the temporal evolution of porosity. In the simplest formulation, α ϕ = 1 , but it is retained here in general form to allow future calibration or nondimensional reformulation.

For a point or line injection source, the well may be represented by an internal singular source term, q ( x , t ) = Q ( t ) δ ( x − x w ) , (5) where Q ( t ) is the injection rate and δ ( · ) is the Dirac delta distribution centered at the well location x w . This representation is particularly convenient for the BEM framework because it avoids volumetric meshing around the well.

The Darcy velocity is governed by v = − k μ ( ∇ p − ρ f g ) , (6)

where k ( x , t ) is the intrinsic permeability, μ is the dynamic viscosity, and 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 v = − k μ ∇ p . (7)

Substituting Darcy’s law into the mass balance equation gives the governing hydraulic diffusion equation S ∂ p ∂ t + α ϕ ∂ ϕ ∂ t − ∇ · ( k μ ∇ p ) = q . (8)

This equation is nonlinear because both ϕ 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 k = k ( x , t ) , x ∈ Ω i , i = 1 , 2 , … , N r , (9) where Ω i denotes the i -th material subdomain. The continuity conditions imposed across internal interfaces are discussed later in this section.

The geochemical alteration of the reservoir is represented through a scalar reaction progress variable ξ ( x , t ) , where ξ = 0 corresponds to the undegraded material state and increasing ξ indicates progressive mineral dissolution. At the Darcy scale, the kinetics may be expressed in a reduced phenomenological form as ∂ ξ ∂ t = R ( p , ξ , χ ) , (10) where R ( · ) is the effective dissolution rate and χ denotes additional thermodynamic or compositional variables, such as CO ₂ concentration, acidity, or saturation state. For the present formulation, and in accordance with the pressure-driven reactive framework adopted in the manuscript, a pressure-dependent first-order kinetic law is introduced: ∂ ξ ∂ t = λ Ψ ( p ) ( 1 − ξ ) , (11) where λ is the kinetic rate coefficient and Ψ ( p ) is a pressure-activation function describing the coupling between local pressure perturbation and reactive driving force. A convenient and numerically robust choice is Ψ ( p ) = p − p 0 p ref , (12) for p ≥ p 0 , with Ψ ( p ) = 0 otherwise , where p ref is a reference pressure scale. Under this definition, the reaction progress equation becomes ∂ ξ ∂ t = λ ( p − p 0 p ref ) ( 1 − ξ ) , p ≥ p 0 . (13)

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 ξ → 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, Ψ ( p ) 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.

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 ∂ ϕ ∂ t = ( ∂ ϕ ∂ ξ ) ∂ ξ ∂ t + ( ∂ ϕ ∂ ε v ) ∂ ε v ∂ t , (14) where ε v = ∇ · u is the volumetric strain.

For the present carbonate dissolution model, the dominant chemical contribution is represented as ϕ = ϕ 0 + a ϕ ξ , (15) or equivalently, ∂ ϕ ∂ t = a ϕ ∂ ξ ∂ t , (16) where ϕ 0 is the initial porosity and a ϕ > 0 is a porosity-increase coefficient. If geomechanical compaction or dilation is retained, the constitutive law may be extended to ϕ = ϕ 0 + a ϕ ξ + b ϕ ε v , (17) where b ϕ is a poromechanical coupling parameter. This form allows both chemical pore enlargement and mechanically induced pore-volume variation to contribute to the storage evolution.

Substituting the reaction law into the porosity relation yields ∂ ϕ ∂ t = a ϕ λ Ψ ( p ) ( 1 − ξ ) + b ϕ ∂ ε v ∂ t , (18) which explicitly shows how pressure, dissolution, and deformation jointly affect pore structure.

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 k = k ( ϕ , σ eff , ξ ) , (19) where σ eff denotes the effective stress tensor or an equivalent stress measure.

A convenient and physically meaningful representation is the multiplicative law k ( x , t ) = k 0 ( x ) exp ( β ξ ξ ) exp ( − β σ Δ σ m ′ ) , (20) where k 0 ( x ) is the initial heterogeneous permeability, β ξ is the reaction-sensitivity coefficient, β σ is the stress-sensitivity coefficient, and Δ σ m ′ is the change in mean effective stress. The factor exp ( β ξ ξ ) captures permeability enhancement due to dissolution, while exp ( − β σ Δ σ m ′ ) accounts for permeability reduction under increased compressive effective stress. When pressure rise leads to stress relaxation in the skeleton, this second term may also increase permeability, depending on sign convention.

If a purely chemical model is desired, the permeability law reduces to k ( x , t ) = k 0 ( x ) exp ( β ξ ξ ) , (21)

Alternatively, a Kozeny–Carman-type relation may be used: k k 0 = ( ϕ ϕ 0 ) m ( 1 − ϕ 0 1 − ϕ ) n , (22) where m and n are empirical exponents. However, the exponential form is especially attractive for nonlinear BEM implementation because it is smooth, monotone, and readily linearized within incremental solution schemes.

To incorporate geomechanical feedback, the porous solid is modeled as a linear poroelastic continuum under quasi-static conditions. Neglecting inertia, the equilibrium equation is ∇ · σ + b = 0 , (23) where σ is the total stress tensor and b is the body force vector. The total stress is related to the effective stress through Biot’s effective stress principle: σ = σ ′ − αp I , (24) or equivalently, σ ′ = σ + αp I , (25) where σ ′ is the effective stress tensor, α is Biot’s coefficient, and I is the identity tensor.

For small strains, the strain tensor is ε = 1 2 ( ∇ u + ∇ u T ) , (26) and the constitutive relation for an isotropic elastic skeleton is σ ′ = 2 G ε + λ s tr ( ε ) I , (27) where G is the shear modulus and λ s is Lamé’s first parameter. Combining these expressions yields the poroelastic equilibrium equation in displacement form, ∇ · [ 2 G ε + λ s tr ( ε ) I ] − α ∇ p + b = 0 . (28)

The volumetric strain relevant for porosity evolution is then ε v = ∇ · u . (29)

To couple mechanics back to flow, the mean effective stress is defined as σ m ′ = 1 d tr ( σ ′ ) , (30) and inserted into the permeability law introduced above. Through this mechanism, injection pressure modifies effective stress, effective stress modifies permeability, and permeability subsequently alters the pressure diffusion process.

Combining the preceding relations, the governing reactive–geomechanical flow model may be written as the following coupled system in Ω × ( 0 , T ] :

2.7.1. Pressure diffusion equation S ∂ p ∂ t + α ϕ ∂ ϕ ∂ t − ∇ · ( k ( ξ , σ m ′ , x ) μ ∇ p ) = q . (31)

2.7.2. Reaction progress equation ∂ ξ ∂ t = λ Ψ ( p ) ( 1 − ξ ) . (32)

2.7.3. Porosity evolution equation ϕ = ϕ 0 + a ϕ ξ + b ϕ ε v . (33)

2.7.4. Geomechanical equilibrium equation ∇ · [ 2 Gε + λ s tr ( ε ) I ] − α ∇ p + b = 0 , ε = 1 2 ( ∇ u + ∇ u T ) . (34)

2.7.5. Permeability constitutive law k = k 0 ( x ) exp ( β ξ ξ ) exp ( − β σ Δ σ m ′ ) . (35)

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.

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: p ( x , 0 ) = p 0 ( x ) , ξ ( x , 0 ) = 0 , ϕ ( x , 0 ) = ϕ 0 ( x ) , k ( x , 0 ) = k 0 ( x ) , u ( x , 0 ) = u 0 ( x ) , (36) where n 0 ϕ 0 and k 0 may be spatially heterogeneous.

2.8.2. Hydraulic boundary conditions

The external boundary Γ may be decomposed into prescribed-pressure and prescribed-flux parts, such that Γ = Γ p ∪ Γ q , Γ p ∩ Γ q = ∅ . (37)

On Γ p , Dirichlet conditions are imposed: p = p ¯ on Γ p × ( 0 , T ] . (38)

On Γ q , Neumann conditions are prescribed: − k μ ∇ p · n = q ¯ n on Γ q × ( 0 , T ] , (39) where n denotes the outward unit normal vector.

For an infinite or semi-infinite reservoir idealization, a far-field condition is adopted: p ( x , t ) → p 0 as | x | → ∞ . (40)

2.8.3. Mechanical boundary conditions

Similarly, the mechanical boundary is partitioned into displacement and traction boundaries: Γ = Γ u ∪ Γ t , Γ u ∩ Γ t = ∅ . (41)

On Γ u , u = u ¯ on Γ u × ( 0 , T ] , (42) and on Γ t , σ n = t ¯ on Γ t × ( 0 , T ] . (43)

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

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 Γ ij separating subdomains Ω i and Ω j .

For hydraulic flow, the continuity of pressure and normal flux requires p i = p j on Γ ij , (44) ( k i μ ∇ p i ) · n ij = ( k j μ ∇ p j ) · n ij on Γ ij , (45) where n ij is the normal vector oriented from Ω i to Ω j .

For mechanics, continuity of displacement and traction is imposed: u i = u j on Γ ij , (46) σ i n ij = σ j n ij on Γ ij . (47)

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.

The constitutive relations above render the problem nonlinear, and therefore an incremental solution strategy is adopted. Let t n + 1 = t n + Δ t . Over each time increment, the variables are updated as p n + 1 = p n + Δ p , ξ n + 1 = ξ n + Δ ξ , ϕ n + 1 = ϕ n + Δ ϕ , k n + 1 = k n + Δ k . (48)

The reaction term is linearized as R ( p , ξ ) ≈ R ( p n , ξ n ) + ∂ R ∂ p | n Δ p + ∂ R ∂ ξ | n Δ ξ . (49) and the permeability update is approximated by Δ k ≈ ∂ k ∂ ξ | n Δ ξ + ∂ k ∂ σ m ′ | n Δ σ m ′ . (50)

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.

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 S ∂ p ∂ t + a ϕ ∂ ξ ∂ t − ∇ · ( k 0 ( x ) e β ξ ξ μ ∇ p ) = q , (51) ∂ ξ ∂ t = λ Ψ ( p ) ( 1 − ξ ) , (52) with ϕ = ϕ 0 + a ϕ ξ , k = k 0 ( x ) e β ξ ξ . (53)

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.

At time level t n + 1 = t n + Δ t , the nonlinear pressure equation in Section 2 is linearized by freezing the coefficients at the previous iteration or previous time level. Let p n + 1 = p n + Δ p , ξ n + 1 = ξ n + Δ ξ , ϕ n + 1 = ϕ n + Δ ϕ , k n + 1 = k n + Δ k . (54)

Over the interval [ t n , t n + 1 ] , the hydraulic operator is approximated using effective coefficients S i ( n ) = S ( ϕ i n , ξ i n , σ m , i ′ n ) , k i ( n ) = k ( ξ i n , σ m , i ′ n , x ) , (55) assumed constant within each subregion Ω i during the current increment. The pressure equation in Ω i then takes the form S i ( n ) ∂ p i ∂ t − ∇ · ( k i ( n ) μ ∇ p i ) = q i ( n ) + f i ( n ) , x ∈ Ω i , (56) where q i ( n ) represents injection or extraction sources and f i ( n ) is an equivalent internal source term collecting the effects of porosity evolution, reaction coupling, and coefficient linearization. A convenient definition is f i ( n ) = − α ϕ ∂ ϕ i ∂ t + r i ( n ) , (57) where r i ( n ) includes residual terms arising from the linearization of permeability and storage. This reformulation converts the original nonlinear problem into a sequence of linear transient diffusion equations with known right-hand-side forcing at each increment. This is precisely the step that makes the reactive model amenable to BEM without resorting to volumetric remeshing.

Introducing the hydraulic diffusivity c i ( n ) = k i ( n ) μ S i ( n ) , (58) the pressure equation may be rewritten as ∂ p i ∂ t − c i ( n ) ∇ 2 p i = 1 S i ( n ) ( q i ( n ) + f i ( n ) ) . (59)

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

For each subregion Ω i , the corresponding adjoint fundamental solution G i ( x , t ; ξ , τ ) is defined as the solution of ∂ G i ∂ t − c i ( n ) ∇ 2 G i = δ ( x − ξ ) δ ( t − τ ) , t > τ , (60) subject to vanishing conditions for t < τ . Here , ξ denotes the source or collocation point, and δ ( · ) is the Dirac distribution.

In two-dimensional infinite space, the transient diffusion fundamental solution is G i ( x , t ; ξ , τ ) = H ( t − τ ) 4 π c i ( n ) ( t − τ ) exp [ − r 2 4 c i ( n ) ( t − τ ) ] , (61) where r = ‖ x − ξ ‖ , (62) and H ( · ) is the Heaviside step function. In three dimensions, the corresponding kernel is G i ( x , t ; ξ , τ ) = H ( t − τ ) [ 4 π c i ( n ) ( t − τ ) ] 3 2 exp [ − r 2 4 c i ( n ) ( t − τ ) ] . (63)

The normal derivative required in the boundary integral equation is given by ∂ G i ∂ n = ∇ G i · n , (64) with n denoting the outward unit normal vector on the boundary. The use of these time-dependent Green’s functions is central to the present formulation, because it permits exact representation of transient diffusion within each subregion while confining spatial discretization to the boundary and internal interfaces.

Multiplying the linearized pressure equation in Ω i by the fundamental solution G i , multiplying the adjoint equation for G i by p i , subtracting the two expressions and integrating over Ω i × [ 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 ξ ∈ Ω i ¯ satisfies c ( ξ ) p i ( ξ , t n + 1 ) + ∫ Γ i ∫ t n t n + 1 ∂ G i ∂ n ( x , t ; ξ , t n + 1 ) p i ( x , t ) dτd Γ = ∫ Γ i ∫ t n t n + 1 G i ( x , t ; ξ , t n + 1 ) q n , i ( x , t ) dτd Γ + ∫ Ω i ∫ t n t n + 1 G i q i ( n ) + f i ( n ) S i ( n ) dτd Ω + ∫ Ω i G i ( x , t n ; ξ , t n + 1 ) p i ( x , t n ) d Ω , (65) where the hydraulic flux is defined as q n , i = − k i ( n ) μ ∇ p i · n , (66) and c ( ξ ) is the standard free-term coefficient, equal to 1 for interior points and to the usual solid-angle fraction for boundary collocation points.

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.

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 x w , the source term is written as q i ( n ) ( x , t ) = Q w ( t ) δ ( x − x w ) , (67) and the associated domain integral reduces immediately to ∫ n t n + 1 G i ( x , t ; ξ , t n + 1 ) q i ( n ) c ( n ) dτd Ω = ∫ n t n + 1 Q w ( τ ) c ( n ) G i ( x w , τ ; ξ , t n + 1 ) dτ , (68)

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 ( n ) , 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 ( n ) 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.

The reservoir is partitioned into N r subregions, Ω = ⋃ i = 1 N r Ω i , (69) with each Ω i treated as locally homogeneous during a given time step, but with coefficients updated between increments. This domain decomposition is particularly effective for representing sharp heterogeneity contrasts, such as layered carbonate facies, damaged zones, or reaction-altered near-well regions.

For each subregion Ω i , a separate boundary integral equation is written over its boundary Γ i = Γ i ext ∪ Γ i int , (70)

where Γ i ext is the external boundary segment and Γ int int denotes interfaces shared with adjacent regions. Across each interface Γ ij = ∂ Ω i ∩ ∂ Ω j , continuity of pressure and flux is imposed: p i = p j on Γ ij , (71) q n , i + q n , j = 0 on Γ ij . (72)

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.

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 and q n , i n + 1 the boundary unknowns at the end of the ( n + 1 ) − th time step. The time convolution integrals involving G i and ∂ G i / ∂ n are then evaluated over the current step and, when needed, accumulated from previous steps through a history term. The generic form becomes c ( ξ ) p i n + 1 ( ξ ) + ∫ Γ i H i ( 0 ) ( ξ , x ) p i n + 1 ( x ) d Γ = ∫ Γ i G i ( 0 ) ( ξ , x ) q n , i n + 1 ( x ) d Γ + b i n + 1 ( ξ ) , (73) where H i ( 0 ) and G i ( 0 ) are the influence kernels integrated over the current time step, and b i n + 1 includes all known contributions from source terms, previous-time history, and nonlinear corrections.

More explicitly, b i n + 1 = b i , hist n + 1 + b i , src n + 1 + b i , nl n + 1 , (74) with the three components corresponding respectively to temporal history, analytical well forcing, and equivalent nonlinear domain effects. The use of implicit time integration stabilizes the transient reactive problem and is particularly beneficial when permeability changes become rapid near the injection zone.

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 : p i ( x ) ≈ p i ( e ) , q n , i ( x ) ≈ q n , i ( e ) , x ∈ Γ i ( e ) . (75)

Substituting these approximations into the boundary integral equation and collocating at the element centers yields c m p i , m n + 1 + ∑ e = 1 N i H me ( i ) p i , e n + 1 = ∑ e = 1 N i G me ( i ) q n , i , e n + 1 + b i , m n + 1 , m = 1 , 2 , … , N i , (76) where the influence coefficients are H me ( i ) = ∫ Γ i ( e ) H i ( 0 ) ( ξ m , x ) d Γ , G me ( i ) = ∫ Γ i ( e ) G i ( 0 ) ( ξ m , x ) d Γ , (77)

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.

Collecting all collocation equations for a subregion Ωi\Omega_iΩi , the discrete system may be written compactly as H i ( n + 1 ) p i ( n + 1 ) = G i ( n + 1 ) q i ( n + 1 ) + b i ( n + 1 ) , (78) or, equivalently, A i ( n + 1 ) x i ( n + 1 ) = b i ( n + 1 ) , (79) where x i ( n + 1 ) contains the unknown boundary values after application of the prescribed external boundary conditions. The matrices H i ( n + 1 ) and G i ( n + 1 ) depend on the updated diffusivity c i ( n ) , and therefore change from one time step to the next as the reaction field evolves.

The right-hand-side vector b i ( n + 1 ) contains: 1.

the transient history contribution from p i n ,

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.

The global reservoir system is obtained by assembling all subdomain equations and imposing interface continuity conditions. Let p = [ p 1 p 2 ⋯ p N r ] T , q = [ q 1 q 2 ⋯ q N r ] T . (80)

After merging common interface unknowns and enforcing p i Γ ij = p j Γ ij , q i Γ ij + q j Γ ij = 0 , (81)

the assembled system assumes the block form A ( n + 1 ) X ( n + 1 ) = B ( n + 1 ) , (82) where X ( n + 1 ) is the global vector of unknown boundary pressures and fluxes over all external boundaries and interfaces.

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.

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 ξ n , ϕ n , and k n , compute the effective storage and diffusivity in each subregion and solve A ( n + 1 ) X ( n + 1 ) = B ( n + 1 ) , (83) for the boundary unknowns. Then evaluate the interior pressure p n + 1 where needed.

3.10.2. Step 2: update the reaction progress variable

Using the pressure field obtained in Step 1, update ξ according to the reaction law introduced in Section 2. For example, with an implicit or semi-implicit Euler rule, ξ n + 1 = ξ n + Δ tλ Ψ ( p n + 1 ) ( 1 − ξ n ) . (84)

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 ϕ n + 1 = ϕ 0 + a ϕ ξ n + 1 + b ϕ ε v n + 1 . (85)

or, in the reduced model, ϕ n + 1 = ϕ 0 + a ϕ ξ n + 1 . (86)

The permeability is updated using k n + 1 = k 0 ( x ) exp ( β ξ ξ n + 1 ) exp ( − β σ Δ σ m ′ n + 1 ) , (87)

or, in the reduced reactive-flow model, k n + 1 = k 0 ( x ) exp ( β ξ ξ n + 1 ) . (88)

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 ‖ k ( m + 1 ) − k ( m ) ‖ ‖ k ( m ) ‖ < ε k , ‖ ξ ( m + 1 ) − ξ ( m ) ‖ ‖ ξ ( m ) ‖ < ε ξ , (89) where ε k and ε ξ are prescribed tolerances. In many practical injection simulations, however, a single staggered update per step is sufficient and numerically stable, especially when Δ t is chosen adaptively.

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.

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.

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 S i ( n ) ∂ p i ∂ t − ∇ · ( k i ( n ) μ ∇ p i ) = q i ( n ) − a ϕ ∂ ξ i ∂ t , (90) with ∂ ξ i ∂ t = λ Ψ ( p i ) ( 1 − ξ i ) , k i ( n ) = k 0 , i exp ( β ξ ξ i n ) . (91)

The corresponding BEM system at time level n + 1 is written as H ( n + 1 ) p ( n + 1 ) = G ( n + 1 ) q ( n + 1 ) + b hist ( n + 1 ) + b well ( n + 1 ) + b react ( n + 1 ) . (92)

Here, b 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.

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 Ω 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 S i ( n ) ∂ p i ∂ t − ∇ · ( k i ( n ) μ ∇ p i ) = q i ( n ) + f i ( n ) , (93) where all coefficients are frozen at the current iterate. For this problem, the time-dependent Green’s function used in Section 3 is the exact fundamental solution of the associated transient diffusion operator. The derived boundary integral equation therefore does not introduce any approximation at the continuous level. The only approximations arise later through temporal discretization, boundary interpolation, numerical quadrature, and nonlinear coefficient updating.

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: (i)

exact Green’s-function representation of the linearized subproblem,

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.

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 Δ t decreases.

Let the nonlinear hydraulic operator be written abstractly as L ( p , ξ , ϕ , k ) = 0 , (94) with constitutive dependencies ϕ = ϕ ( ξ , ε v ) , k = k ( ξ , σ m ′ , x ) . (95)

Over a time step, the operator is approximated by L ( p n + 1 , ξ n + 1 , ϕ n + 1 , k n + 1 ) ≈ L ( n ) ( p n + 1 ) + R ( n ) , (96) where L ( n ) denotes the frozen-coefficient linear operator and R ( n ) collects residual terms associated with constitutive evolution. If the constitutive functions are smooth and the update is performed consistently, then R ( n ) = O ( Δ t ) , (97) for a first-order incremental scheme. Hence, as Δ t → 0 , the linearized subproblem converges to the original nonlinear problem in the usual local truncation sense.

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.

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 λ = 0 , then ξ ( x , t ) = 0 , ϕ = ϕ 0 , k = k 0 , (98) and the model reduces to the classical transient diffusion equation S ∂ p ∂ t − ∇ · ( k 0 μ ∇ p ) = q . (99)

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 λ = 0 but k 0 ( x ) varies between subregions, the formulation reduces to transient diffusion in a piecewise heterogeneous reservoir. The exact enforcement of p i = p j , q n , i + q n , j = 0 . (100) across every interface ensures that the method reproduces the correct transmission of pressure waves through permeability contrasts. Agreement with reference solutions in this setting verifies the subdomain assembly procedure and the interface-coupling strategy.

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 k ( t ) = k 0 exp ( β ξ ξ ( t ) ) . (101) may be compared with a semi-analytical ordinary differential reduction of the governing system. Recovery of this limit confirms that the nonlinear constitutive updating has been implemented correctly and does not introduce artificial damping or amplification.

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.

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 ref ( x , t ) denote the reference solution and p BEM ( x , t ) the numerical result. The pointwise relative error may be defined as e p ( x , t ) = | p BEM ( x , t ) − p ref ( x , t ) | max ( | p ref ( x , t ) | , ϵ ) , (102) where ϵ is a small regularization parameter preventing division by zero. Global accuracy measures may then be written as E L 2 ( t ) = ( ∑ m = 1 N o [ p m BEM ( t ) − p m ref ( t ) ] 2 ∑ m = 1 N o [ p m ref ( t ) ] 2 ) 1 2 , (103) and E ∞ ( t ) = max 1 < m < N o | p m BEM ( t ) − p m ref ( t ) | , (104) where N o is the number of observation points.

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.

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 ‖ p − p h ‖ → 0 as N b → ∞ (105) where p h is the BEM approximation. For constant boundary elements, first-order convergence in an appropriate norm is generally expected, while higher-order boundary interpolation yields correspondingly improved rates for smooth solutions.

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

pressure at the wellbore,

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 η h = | J h − J h / 2 | | J h / 2 | . (106)

Convergence is established when η 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 , ξ , 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.

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 ‖ p Δ t − p ‖ = O ( Δ t ) , ‖ ξ Δ t − ξ ‖ = O ( Δ t ) , (107) provided that the constitutive updates remain sufficiently smooth.

A temporal refinement study is performed by computing the solution with decreasing time steps Δ t , Δ t 2 , Δ 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: p ( x , t ) remains finite for finite t , (108) 0 ≤ ξ ( x , t ) ≤ 1 (109) ϕ ( x , t ) > 0 , k ( x , t ) > 0 . (110)

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.

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 Γ ij , the pressure mismatch may be quantified by E p Γ ij = ‖ p i − p j ‖ L 2 ( Γ ij ) max ( ‖ p i ‖ L 2 ( Γ ij ) , ϵ ) , (111) while the flux imbalance may be measured by E q Γ ij = ‖ q n , i + q n , j ‖ L 2 ( Γ ij ) max ( ‖ q n , i ‖ L 2 ( Γ ij ) , ϵ ) . (112)

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 ∫ Ω ∂ ∂ t ( ϕ ρ f ) d Ω + ∫ Γ ρ f v · n d Γ − ∫ Ω q m d Ω ≈ 0 . (113)

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.

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: 1.

a known analytical radial solution in a homogeneous medium,

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.

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 k n + 1 = k 0 exp ( β ξ ξ n + 1 ) . (114) or, in the more general case, k n + 1 = k 0 exp ( β ξ ξ n + 1 ) exp ( − β σ Δ σ m ′ n + 1 ) . (115)

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 E k n + 1 = ‖ k n + 1 − k n ‖ ‖ k n ‖ , (116) together with the residual between successive internal nonlinear iterations, R k ( m + 1 ) = ‖ k ( m + 1 ) − k ( m ) ‖ ‖ k ( m ) ‖ . (117)

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.

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: 1.

pressure histories at monitoring points,

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.

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.

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.

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

recovery of classical transient diffusion solutions in nonreactive limits,

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.

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.