The effects of spatially correlated heterogeneities on acidising wormhole development

A_v

dimensionless interfacial area per unit volume available for a reaction

a₀

initial interfacial area per unit volume available for a surface reaction

a_v

interfacial area per unit volume available for a surface reaction

C₀

inlet concentration of the acid

C_f

acid concentration in bulk fluid

C_s

acid concentration at the solid–fluid interface

$Cφ0φ0$

anisotropic covariance matrix

c_f

dimensionless acid concentration in bulk fluid

c_s

dimensionless acid concentration at the solid–fluid interface

Da

Damköhler number

D_e

effective dispersion tensor

D_m

molecular diffusion coefficient

H

dimensionless effective separation distance

$HT2$

macroscale Thiele modulus

h_i

component in the direction i of the effective separation vector

$hT2$

pore-scale Thiele modulus

K

permeability tensor

K₀

initial average permeability of rock

k_c

mass transfer coefficient

k_s

surface reaction rate constant

L

length of the core

N_ac

dimensionless acid capacity number

P

pressure

P_out

pressure at the exit boundary of the core

Pe

Péclet number

PV_BT

pore volume required to break through

p

dimensionless pressure

R_s

surface reaction rate

R_t

transport rate by mass transfer (diffusion) from the bulk to the surface–fluid interface

r

dimensionless mean pore radius

r₀

initial mean pore radius

r_p

mean pore radius

Sh

Sherwood number

t

time

t*

dimensionless time

U

velocity

u

dimensionless velocity

u₀

inlet (injection) velocity

u_t

dimensionless total Darcy flux

x

dimensionless axial/flow direction

x*

axial/flow direction

y, z

transverse directions

y*, z*

dimensionless transverse directions

α

dissolving power of acid

α_0s

coefficient of molecular diffusion in diffusion dominant regimes

β

pore broadening parameter

Δp₀

initial pressure difference across the core

Δp_t

pressure difference across the core at time t

Δφ₀

random porosity heterogeneity range

η

ratio of mean pore diameter to core length

κ

magnitude of dimensionless permeability

λ_i

correlation length scales in direction i

$μφ0$

mean of the correlated φ ₀ field

ξ_T

coefficient of transverse hydraulic dispersion in advection-controlled regimes

ξ_X

coefficient of longitudinal hydraulic dispersion in advection-controlled regimes

ρ_s

solid phase density

$σφ02$

variance of φ ₀

φ

evolving porosity of the rock

φ₀

initial mean porosity

Rock matrix acidisation is a process by which the sediments, mud solids and minerals within the pores of siliciclastic or carbonate reservoirs are dissolved by injecting an acidic solution into the rock below the fracturing pressure level. The rock-matrix-acidising process can enhance the injectivity and productivity of fluid and gas recovery operations in applications including hydrocarbon or geothermal enhanced recovery. The process features a strongly coupled flow and geochemistry, where the physical properties of rock such as porosity and permeability dynamically evolve across multiple time and length scales.

Reactions between chemical solutions and the solid phase appear in various engineering and natural subsurface processes, including irrigating water discharge to aquifer systems (Valdes-Abellan et al., 2017), karst formation (Evans and Lizarralde, 2003; Zhao et al., 2013), chemical interaction in saturated and unsaturated soils (Atchley et al., 2014; Chen et al., 2009; Cubillas et al., 2005; Ouhadi et al., 2006), radioactive waste containment and remediation (Sedighi et al., 2015, 2018; Spycher et al., 2003), geologic sequestration of carbon dioxide (e.g. (Islam et al., 2014, 2016; Kampman et al., 2014; Rochelle et al., 2004)), acid injection for enhanced oil recovery (e.g. the paper of Ghommem et al. (2015) and references therein) and acid-based enhanced geothermal system (Na et al., 2016; Portier and Vuataz, 2010; Portier et al., 2007; Xiong et al., 2013).

Over a span of decades, a compelling body of works has focused on developing transport modelling tools, tested/validated against experiments to predict wormhole formation, propagation and stability in the context of oil recovery (Bastami and Pourafshary, 2016; Cohen et al., 2008; Daccord et al., 1989; Fredd and Fogler, 1998; Ghommem et al., 2015; Glasbergen et al., 2009; Golfier et al., 2006; Kalia and Balakotaiah, 2007; Liu and Liu, 2016; Liu et al., 2017; Maheshwari et al., 2013; Nierode et al., 1972; Qiu et al., 2013; Schechter and Gidley, 1969). In petroleum engineering, acid solutions (typically 15 wt% (Hung et al., 1989)) are used to restore or enhance oil production that is affected by the particles that migrated along the flow and eventually accumulated near the wellbore region during the production lifetime. The preferential dissolution pathways are generated only if certain conditions of rock, acid and injection exist. The parameters and conditions of rock that control efficiency include mineralogical configuration and distribution and petrophysical properties and distributions of rock. The right acid type, volume, dissolving power and injection rate are also controlling properties for efficient wormhole generation. Therefore, the interplay between flow and geochemistry determines the states at which wormholes are generated. The interplay is analysed by describing the parameters of flow and geochemistry through the Péclet number (Pe) and Damköhler number (Da).

Porosity and permeability heterogeneity in the carbonate reservoir play an important role in the formation of wormholes during the reactive dissolution process of rock acidisation. The naturally occurring heterogeneities lead to an uneven permeability increase at the solid–liquid interface. The heterogeneity of reservoir porosity and permeability has been commonly described with a random number in previous works that have dealt with modelling of wormhole formations (see e.g. the papers of Panga et al. (2005), Kalia and Balakotaiah (2009) and Maheshwari and Balakotaiah (2013a)). There are very limited studies that have implemented a realistic distribution of petrophysical properties. The most important difference between a realistic distribution of petrophysical properties and a random distribution is the presence of a spatial correlation with different length scales in the porous system. Liu et al. (2012) adopted a normal distribution within a two-dimensional (2D) radial system with circumferential and radial correlations. A spatial correlation through sequential Gaussian simulation has been used by Liu et al. (2016) to generate the initial porosity distribution. In neither of these works has the influence of varying/increasing correlation lengths in the model been investigated – that is, if a comparison has been made between uniform distribution porosity and permeability, a comparison between varying correlation lengths has not been made. Moreover, the impact of correlation length on the development of wormholes and in particular the optimal pore volume (PV_BT) is not clarified by Liu et al. (2016) (see Figure 16 of their paper).

The study presented here aims to address the aforementioned shortcoming in literature – that is, how different values of correlation length impact the wormhole development. The authors specifically look at the simplified reaction of calcite with hydrochloric acid (HCl) for a fully homogeneous calcite system, but with heterogeneous rock properties (porosity/permeability). All the assumptions made throughout this paper were justified previously by Panga et al. (2005). The authors use an in-house simulation code (based on the two-scale model of Panga et al. (2005)). In the following, a critical literature review of the modelling approaches using random or correlated porosity and/or permeability is first described (see the section headed ‘Heterogeneities of rock and modelling approaches’). Then, the development of a reactive transport model and the governing equations used in the simulation of rock matrix acidisation in this study and the algorithm for generating correlated random fields are described (see the section headed ‘Methodology and the simulation model’). The results are presented and discussed in the section headed ‘Results and discussions’, followed by the section headed ‘Conclusions and future works’.

Early-stage simulations of porosity heterogeneity have been performed using the average continuum model for carbonate dissolution which is based on the Darcy scale and the pore scale for a linear flow (Panga et al., 2005). The initial porosity field consists of generation of random values fluctuating between a range around the mean value of porosity, which are added to the average porosity – that is, [φ₀ − Δφ₀,φ₀ + Δφ₀]. This method is referred to as uniform distribution method that considers a magnitude of randomly generated heterogeneity (Δφ₀/φ₀) within the domain. Permeability is directly related to porosity using semi-empirical relations (e.g. cubic relationship). Kalia and Balakotaiah (2007) extended this modelling approach for a radial flow and showed that a critical value of randomly generated heterogeneous porosity exists that below which the minimum pore-volume-to-breakthrough (PV_BT) is much higher than that of a more heterogeneous porous medium. The critical value indicates the minimum amount of heterogeneity required to create instability in the dissolution front and transition from the planar face dissolution to the formation of wormholes. Kalia and Balakotaiah (2009) conducted an investigation on the effect of medium heterogeneities, such as magnitude and length scale over which the porosity varies, and their effect on reactive dissolution for linear flow in a three-dimensional (3D) domain. They showed that heterogeneity affects the structure of the dissolution patterns and, more importantly from a practical point of view, the amount of acid required to achieve a given increase in permeability. They suggested that an optimum heterogeneity magnitude value at which the least amount of PV_BT is required may exist. Additionally, they demonstrated that the amount of acid required to break through is dependent on the initial rock porosity and dimensions of rock being acidised.

Cohen et al. (2007) proposed a dual-porosity approach using the average-continuum model and considered the reactive medium to consist of two media with different porosities. One medium contains the dominant growing wormholes, and within the other medium, only the short-lived wormholes (and compact dissolution) can occur. Cohen et al. (2007) presented the first 3D simulations of the radial flow, and their model is able to simulate the dissolution patterns known from experiments. Both Cohen et al. (2008) and Kalia and Balakotaiah (2009) showed that the optimum rate or minimum PV_BT was not affected by the domain size.

Liu et al. (2012) introduced the heterogeneity of porosity to be distributed normally as opposed to be distributed uniformly. In the 2D average continuum radial flow model presented by Liu et al. (2012), it is shown that when considering normally distributed porosities, less PV_BT is required compared with when the uniform distribution method is applied. In addition, the optimum injection velocity was found to be lower. Acid flows into high-porosity regions, causing wormholes to develop more unevenly compared to the uniform distribution method. It is shown that a large increase in the circumferential correlation length leads to thicker wormhole geometry and greater diversion of propagation direction of the wormholes. A large increase in the radial correlation length causes formation of an increased number of wormholes and thinner wormholes, which may actually reduce the efficiency of the operation. A critical value for the standard deviation (i.e. range of porosities) was derived by Liu et al. (2012). Above this value, PV_BT decreases sharply, and below this value, it is insensitive. PV_BT is least when the standard deviation is unity.

The average continuum modelling approach was extended to 3D simulations by Cohen et al. (2008), Maheshwari and Balakotaiah (2013b) and De Oliveira et al. (2012). The effect of large-scale porosity–permeability heterogeneity has been investigated through assessment of vuggy-type carbonates (where cave-like porous features exist in the porous media) and fractures (Izgec et al., 2010; Kalia and Balakotaiah, 2009; Liu et al., 2012). They showed that when incorporating vugs, the acid propagates through the wormholes at a much faster rate than in homogeneous rocks. Acid follows a preferential flow path, guided by the vug network, leading to a decrease in PV_BT. Their model, however, does not consider the acid–rock reaction mechanism (Izgec et al., 2008). Later, Izgec et al. (2010) demonstrated through using Darcy–Brinkman formulations for flow that PV_BT is influenced by the connectivity of vuggy pores and its spatial distribution in the domain. Maheshwari and Balakotaiah (2013a) demonstrated that for the 3D simulations, a critical porosity heterogeneity (Δφ₀) exists at which PV_BT is minimum. An increase in the branching and fractal nature of the wormholes is apparent when the randomly generated heterogeneity magnitude (Δφ₀) is increased (Maheshwari and Balakotaiah, 2013a).

More recently, Liu et al. (2016) investigated the effect of incorporating the correlation spatial distribution of the rock petrophysical properties. Lithologically, the authors considered a uniformly distributed porosity of 0·2 ± Δφ₀ and three cases of correlated porosity fields with average porosities of 0·2038, 0·2010 and 0·2029 and similar unreported correlation lengths. The authors considered two different rock-types (dolomite and limestone), with different solubility but with similar porosity/permeability and that are located close to each other. Mineral heterogeneity is established through simulation of grid cells with different solid volumes. The initial porosity distribution is spatially correlated through the use of sequential Gaussian simulation. The authors demonstrated that the location of wormhole propagation paths is highly influenced by the locality of spatially correlated high-porosity features. However, as stated in the section headed ‘Introduction’, the correlation length is kept constant (for all the correlated cases) and, in practice, the comparison is based on one random distribution and three correlated cases with a similar spatial correlation length. Moreover, the conclusion made in the paper that ‘that larger pore scale heterogeneity leads to a smaller PV_BT’ is not well established in results (Liu et al., 2016: p. 91; Figure 16). Therefore, using variably correlated heterogeneities, the authors investigate how different correlation lengths may affect wormhole development and PV_BT, in an approach more systematic than that of Liu et al. (2016).

Based on the two-scale average continuum model and a dimensionless form of the governing equations, a two-phase reactive transport model has been developed (Babaei and Sedighi, 2018). The formulation facilitates a comparative setting for obtaining the ideal Pe–Da regime at which the effective wormholing can be achieved (see Figure 1). The governing equations are presented briefly here. Details of the theoretical development and numerical model can be found in the paper of Babaei and Sedighi (2018). In deriving these equations, most importantly the reaction of calcite with HCL solution is simplified into experimentally driven linear kinetics (Alkattan et al., 1998). These simplifications have previously been used successfully for carbonate matrix acidisation modelling validated by experimental results (Ghommem et al., 2015).

where R_s is the surface reaction rate which is a function of the reaction rate constant, activity of H⁺ (⁠$γH+,s$⁠) and surface concentration of H⁺ (⁠$CH+,s$⁠), and R_t is the transport rate by mass transfer (diffusion) from bulk to the surface–fluid interface. At the stationary state, the quantity of hydrogen ions transported to the surface equals that consumed by the reaction (Alkattan et al., 1998). Assuming that the exponent of the reaction rate with respect to the concentration of acid is ∼1

where k_ck_s/(k_s + k_c) is the overall (observed) calcite dissolution rate constant that is independent of solution composition, but to a large extent dependent on the experimental design (Alkattan et al., 1998). Therefore, there are uncertainties in this variable that are embedded into the range of Damköhler numbers in this study.

To use the dimensionless equations, a series of dimensionless variables are defined as reported in Table 1. Despite the variations in the mass transfer coefficient, the authors use a constant value (Sh_∞ = 3). Therefore, in the dimensionless analysis, the Sherwood number reflects the value of the mass transfer coefficient. Similarly, the Damköhler number is a reflection of k_s.

Based on the two-scale formulation of Panga et al. (2005), the dimensionless flow and transport equations can be simplified to

The dimensionless porosity development equation is written as

The following semi-empirical equations are used between the porosity, permeability, average pore radius and surface area of the rock

These relationships were used by Panga et al. (2005) and Maheshwari et al. (2013) to model carbonate rocks acidising. If β = 1, the permeability evolution reduces to the well-known Kozeny–Carman correlation by Wyllie and Gardner (1958): $κ∝φ3/(1−φ)2$⁠. Panga et al. (2005) extended the Kozeny–Carman correlation to a dissolving medium by including β in the relationships.

The correlated porosity–permeability fields are generated using the field generator developed by Nowak et al. (2008). The authors used fast Fourier transform based on the power spectral estimation method that estimates the spectral density function from a random autocorrelated field. The anisotropic covariance matrix for a second-order stationary field is defined using exponential, Gaussian or spherical methods (Nowak et al., 2008). In this work, the authors use the exponential covariance matrix to produce correlated maps of φ₀, given as

where $σφ02$ is the variance of φ₀ and H is the dimensionless/effective (anisotropic) separation distance scaled by the correlation length scales λ_i[L], i = x, y, z: $H=(hx2/λx2+hy2/λy2+hz2/λz2)0⋅5$⁠, where h_i[L] is the separation vector component in the direction i.

The mean of φ₀, $μφ0$⁠, is added to the fast Fourier transform of the random autocorrelated field generated from $Cφ0φ0$ that has a zero mean.

The variables required for modelling the wormhole development based on the dimensionless equations are assigned values as reported in Table 2. The values are compared with experimental data from core experiments of Ghommem et al. (2015) for 15% HCL solution injected into carbonate rock at 65°C and with 1 ml/min < u₀ < 25 ml/min.

The HCL with dissolving power of N_ac is injected into a purely homogeneous calcite mineral which its dissolution by acid is governed by Da everywhere. The three Da values of 10, 100 and 1000 are used to represent change in the reaction rate constant (k_s) or initial interfacial area per unit volume available for surface reaction (a₀). This range of values for Da agrees with the experimental data provided by Ghommem et al. (2015) (see Figure 5 of their paper). There, the authors change u₀ to have 0·1 < Da < 10⁶.

The correlated porosity fields are generated using $μφ0$ = 0·27 and $σφ02$ of 0·0045, so that the mean of porosity is comparable with those of existing simulation studies (Panga et al. (2005), Ghommem et al. (2015) and Liu et al. (2016): $μφ0$ = 0·20). Two correlation lengths (λ_x,y,z) of two and five gridblocks (i.e. 0·0667 and 0·1667 cm) are used to generate the correlated fields. The correlation lengths in each direction are equal – that is, λ_x = λ_y = λ_z – so that the field is isotropically correlated. To make a comparison, the authors also generate a random distribution of porosity with $μφ0$ of 0·27 and $σφ02$ of 0·0045. Therefore, in total, three porosity fields are considered for comparison as shown in Figure 2. The authors note that, in this figure, only porosities larger than 0·3 are shown to make the correlated features inside the domain visible. Also, the authors note that the axes’ labels refer to the gridblock numbers and not the physical lengths scales.

The histograms of the porosity fields under comparison are shown in Figure 3. There is a clear difference between the distribution of the uncorrelated random field and correlated fields. The correlated fields demonstrate a realistic representation of porosity (e.g. Figure 5 of the paper by Sahin et al. (2003)) for the porosity distribution for an Upper Jurassic carbonate reservoir located in the Eastern Province of Saudi Arabia). As such, the uniform distribution of a random porosity field – which is commonly used in the literature of acidisation simulation – is not realistic. The three initial porosity fields make the pore volume (PV) of three simulations equal to 5·003, 5·5479 and 5·8085 cm³. The difference between these values will not make the authors’ comparisons inconsistent, because they use pore volume injected to create a non-dimensional injection rate.

Using the porosity fields, the authors generate initial permeability fields (K₀) using a simple scaling so that the permeabilities are between 0 and 100 mD (9·8692 33×10⁻¹⁴ m²). A constant inlet velocity of u₀ is used which is equal to a constant flow rate of 1 cm³/s divided by the cross-sectional area of the core sample, which is 4 cm². Therefore, u₀ is 0·25 cm/s. Using an effective molecular diffusivity of acid equal to 3 × 10⁻⁹ m²/s (Liu et al., 2016), the macroscopic Péclet number is 4·17 × 10⁴ in the authors’ simulations. The large Péclet number lies on the upper limits of experiments carried out by Panga et al. (2005); therefore, the generation of wormholes will depend only on the magnitude of the Damköhler number.

The finite-difference and implicit pressure and explicit saturation schemes are employed to solve the pressure and continuity equations (Equations 4 and 5) and the concentration equation (Equation 6) (Babaei and Sedighi, 2018). A 150 × 60 × 60 domain (gridblock dimensions are 0·0333³ cm³) is used to discretise the governing equations. The criterion for the conditions in which wormholes can be considered to have been developed fully is not universally defined in different literature studies. While many authors use a greater-than-a-threshold breakthrough concentration of acid as an indication of wormhole development, one can argue that this may happen in non-wormholing dissolution regimes as well. To this end, the authors use the pressure gradient across the domain in the flow direction as a practical and reliable criterion. More specifically, they use the ratio of pressure difference at any simulation time t over the initial pressure gradient (Δp_t/Δp₀) as a measure of how effectively wormholes could reduce the pressure difference. If the dimensionless pressure drop (Δp_t/Δp₀) has become less than or equal to 10⁻⁴, the authors assume the wormholes are generated and have fully grown to create superconductive flow paths so that the pressure difference has dropped to almost zero across the domain. The simulations are continued as long as either the dimensionless pressure drop reaches the threshold of 1 × 10⁻⁴ or the pore volume of injected acid (PVI) reaches 4. Both these values are chosen arbitrarily/heuristically, but as will be shown in the next section, these values lead to effective identification of wormhole generation.

The authors base their comparisons on nine simulation scenarios: three porosity–permeability fields and three Damköhler regimes of Da = 10, Da = 100 and Da = 1000. In all nine simulation cases, thanks to a high Péclet number, wormholes are generated (no face or uniform dissolution patterns are observed). Figure 4 shows how wormholes are generated for the simulation scenarios. There is clearly a pattern shift from a low Da of 10 to a Da of 1000. The wormholes are ramified for a low Da, and the porosity increase within the wormholes is around 0·5. Increasing Da to 100, the wormholes are sharpened (so that a reduced amount of acid is consumed), yet still ramified. Finally, for Da = 1000, the wormholes are sharp and effective, consuming the least amount of acid. In terms of comparing the wormholes between the different porosity fields, the authors have to first note that Figures 4(a)–4(c) are produced at PVI = 4 (because Δp_t/Δp₀ does not reach 10⁻⁴), whereas Figures 4(d)–4(i) are produced exactly when the dimensionless pressure drop has reached 10⁻⁴. Therefore, the profiles shown in Figures 4(d)–4(i) are not generated for the same amount of acid injected and show only how the wormholes have developed when the simulations are actually terminated.

In Figure 5, the amount of acid consumed in terms of PVI is plotted for different simulation scenarios. The figure clearly shows that moving towards the correlated porosity fields, the PVI required to hit the threshold of Δp_t/Δp₀ = 10⁻⁴ is decreased for Da = 100 and Da = 1000. In Table 3, decreases in PVI for Da = 100 and Da = 1000 are recorded. One can observe that the correlated field with λ_x,y,z = 0·1667 cm has decreased the PVI by 12·2% for Da = 100 and by 17·6% for Da = 1000 relative to the random field. These relative reductions are more considerable compared to the relative reductions in PVI when the randomly generated heterogeneity (Δφ₀/Δφ₀) is increased from 0·5 to 1·5 (see Figure 6 of the paper by Panga et al. (2005)). There, the authors showed that the increased heterogeneity will have a maximum 10% reduction in the PVI for a wormhole regime of Da = 1000 (see Figure 7 of the paper by Panga et al. (2005)).

The reason why the correlated fields facilitate a faster and more effective growth of wormholes is that the correlated porosity–permeability allows either fewer branches to develop from the wormholes or fewer unsuccessful wormholes to generate. The spatial correlation in porosity and permeability favours acid front movement longitudinally and suppresses to some extent its lateral movement (branching). This is more obvious when comparing the porosity increase snapshots of Figure 4(d) with Figure 4(f) for Da = 100 or Figure 4(g) with Figure 4(i) for Da = 1000.

To inspect the effect of correlation length on wormhole geometry (branchiness, thickness and straightness), in Figure 6, the ratio of permeability (κ) is compared between the random porosity field and the correlated porosity with λ_x,y,z = 0·1667 cm, for Da = 100 and Da = 1000 when Δp_t/Δp₀ = 10⁻⁴. Similar to the porosity increase, the snapshots show that the following.

• The sizes of branches for the random porosity field are slightly larger than the correlated porosity field for both Da = 100 and Da = 1000.

• The number of undeveloped wormholes for the random porosity field are slightly larger than the correlated porosity field for both Da = 100 and Da = 1000.

Therefore, both variables – that is, the porosity increase and the ratio of dimensionless permeability over initial permeability – indicate that for a small system of 5 cm, the effects of correlation length on petrophysical properties are tangible. It should be noted that the effect of correlation in the porosity/permeability and the relative reduction in PV_BT (or PVI to make Δp_t/Δp₀ = 10⁻⁴) should be higher for larger model samples – for example, when L = 50 cm.

The findings of this work have clear practical relevance to control and design of acidising operations where the aim is to enhance hydrocarbon or heat recovery from hydrocarbon or geothermal reservoirs. An inevitable existence of correlated features in porous media of these systems means that an acid operation can be carried out more efficiently by consuming less acid. For example, for a target area around the wellbore, 17·6% less PVI translates into several cubic metres less acid consumed. Moreover, ignoring the spatial correlation may result in overconsuming and jeopardising the integrity of reservoir barriers or cap rock, consequently leading to potential contamination of underground water resources by acidic solutions. So far, the objective of most of the acid wormholing modelling works has matched the core sample experiments. However, the spatial correlation of porous media petrophysical properties – that span beyond the core and mesoscopic scale into reservoir scale – should be incorporated in modelling-based operation design.

The acid treatment studied in this work did not include the effects of secondary phases/fluids remaining from a pre-flush stage. However, recent publications, including that of Babaei and Sedighi (2018), have investigated the presence of such phases (e.g. oil blobs or trapped air) on the wormholing process. The investigation (on uniformly distributed porosity fields) shows that the secondary phase will lead to less branchy wormholes and as a result will enhance wormhole development − the results that are supported by experiments (e.g. Shukla et al., 2006). An investigation to depict the interplay of realistic heterogeneity of porous rock (with spatial correlation) and two-phase flow conditions is a of matter of interest in future studies.

Finally, the mechanical stability of developed wormholes and their dependence on the correlation lengths of the porous system are very important in the design of the acidising operation. The process of effective wormhole generation and its stability during the production period are a coupled phenomenon that requires studying the reactive transport and mechanics of the problem alongside each other. During the generation of wormholes, pore fluid pressure withstands the overburden pressure. However, excessive lateral growth of the wormholes and loss of acid solution to rock could generate geomechanical instability. Under this condition, the wormholes become susceptible to collapse. In this line of research, optimising acid reaction rates, amount, acidity and level of interaction with rock are key factors in obtaining the desired effects on the formation at downhole conditions. Coupling of thermohydromechanical effects (e.g. Salimzadeh et al., 2018) is required to determine optimal condition for sustainable acidisation.

The development of the 3D continuum scale code for reactive transport processes was achieved by the support of Engineering and Physical Sciences Research Council First Grant EP/R009678/1 that was awarded to MB.