Simulation and Prediction of Injection Pressure in CO 2 Geological Sequestration Based on Improved LSO Algorithm and BP Neural Networks

In order to avoid the occurrence of caprock integrity damage and gas escape due to injection pressure overrun in CO 2 sequestration, an optimized back propagation (BP) neural network model based on Monte Carlo simulation (MCS) and an improved lion swarm optimization (ILSO) algorithm is proposed for the maximum sustainable injection pressure prediction. Firstly, a hydromechanical model is constructed to simulate the damage changes of the reservoir caprock during injection by ABAQUS. Secondly, in view of the uncertainties of formation parameters that could lead to deviations between the model calculation results and actual geological condition, the MCS method is used, and then, the probability distribution interval of the maximum injection pressure with high probability of caprock failure under di ﬀ erent formation parameters is obtained by MATLAB. This e ﬀ ectively reduces the uncertainty and improves the calculation accuracy. Finally, based on the numerical simulation results, the maximum injection pressure prediction model is constructed. Aiming at the problems of the sensitivity of the BP neural network to initial weights and its poor convergence, tent chaotic mapping and di ﬀ erence mechanism are introduced to improve the LSO algorithm. Following this, the neural network is optimized by ILSO algorithm whose superiority is veri ﬁ ed through 8 benchmark functions, and the maximum injection pressure is e ﬀ ectively predicted according to porosity, permeability, and other parameters. Experimental results show that, compared to the other three optimization methods, the ILSO_BP model has a faster convergence speed, higher prediction accuracy, and stability, which can provide a powerful guide for the safe injection of CO 2 and e ﬃ cient sequestration management.


Introduction
High CO 2 concentration in the atmosphere can lead to a series of problems such as global warming and extreme climate, thus posing a serious threat to the sustainable development of the ecological environment. Geological CO 2 sequestration is one of the most scientific and effective means for controlling the CO 2 concentration in the atmosphere, which is of great significance to mitigate the green-house effect, improve oil and gas recovery, and promote the development of low-carbon economies [1]. The injection of large amounts of CO 2 into the geological reservoir will break the equilibrium in the reservoir and have an important impact on the mechanical stability of the caprock [2] ( Figure 1). Once a leak occurs, it will likely cause a serious environmental disaster. Thus, it is necessary to quantitatively characterize the extent of the caprock damage during the injection. Many complex factors affect the caprock integ-rity in CO 2 sequestration projects, such as the injection rate and material properties [3][4][5]. Among them, the injection pressure is the key external factor triggering the geomechanical effect, whereas the mechanical properties of rocks are the objective internal factor affecting the sealing performance. Moreover, the injection pressure directly affects the efficiency, storage capacity, and potential hydraulic fracture risk. Therefore, exploring an efficient and accurate prediction and evaluation method for calculating the maximum sustainable injection pressure (hereafter referred to as the maximum injection pressure) is a particularly important basis for ensuring the safety of the CO 2 sequestration project, improving the storage efficiency, and reducing the cost.
Several scholars have carried out studies on the maximum injection pressure in CO 2 sequestration, including theoretical analysis methods [6,7], laboratory tests or gas injection tests in shallow sites [8], and numerical simulations [9]. Among them, the numerical simulation method can efficiently make up for the shortcomings of other methods ( Figure 2) and has been widely used in many aspects such as potential assessment of CO 2 sequestration sites [1], analysis of the formation pressure change during injection [10], leakage risk assessment [11], and analysis of the influence of injection parameters in CO 2 enhanced oil recovery [12]. Geological CO 2 sequestration is a complex thermal-hydrological-mechanical-chemical coupling process, and the constitutive models which are used for simulating this process are mainly of three types: elastic constitutive models, plastic constitutive models, and discontinuous media models [13]. The impact of different constitutive models on the simulation results of the caprock integrity analysis is different. Related researches show that the key to ensuring the integrity of the caprock is to make sure that the reservoir pressure does not exceed the breakthrough pressure of the caprock [14], whereas in engineering application, they often increase the injection pressure to improve the CO 2 storage efficiency. The maximum safe injection pressure in the traditional method is obtained by multiplying the formation breakthrough pressure by a safety factor (0.8-0.9). However, this method can only roughly determine a range and cannot accurately predict whether the injection pressure will exceed the critical value and lead to the risk of rupture [15,16]. In addition, the mechanical and permeability properties of reservoir caprock are obtained from the core test data, and porosity and permeability, as important parameters affecting the sequestration performance, present a complex multiparameter, nonlinear mapping relationship with the maximum injection pressure [17]. However, the measurement and sampling errors and the complexity of the geological environment lead to the incompleteness and some uncertainty in the logging data [18,19]. The mean value of the formation parameters in the field is generally used in the model calculation, which leads to the deviations between the calculated results and the actual geological conditions. This, to some extent, reduces the reliability of the model and increases the risk of the sequestration project. Therefore, it is necessary to perform an uncertainty analysis on the relevant formation parameters to accurately predict the maximum injection pressure. However, at present, there are relatively few studies on this aspect.
Artificial neural networks (ANNs) have unique advantages in the analysis and processing of massive, highdimensional big data and can efficiently characterize the complex nonlinear mapping relationships between the input and output data. As a result, ANNs have been widely used in mining the potential laws and predicting the unknown data in the field of oil and gas development, e.g., the prediction of the fracturing capacity of tight sandstone reservoirs based on the genetic algorithm-Elman neural network [20], the remaining oil plane distribution prediction based on machine learning [21], and the seismic reservoir prediction based on a convolutional neural network [22]. In particular, ANNs have also been widely used in CO 2 sequestration engineering, e.g., deep learning-based reservoir mechanical response at the storage site [23], predictive characterization  2 Geofluids of the subsurface porosity by fusing ANNs [24], and ANNbased trap effectiveness prediction in CO 2 sequestration in saline aquifers [25]. Related studies have achieved better results in the fine characterization of reservoirs, evaluation of CO 2 storage safety, and prediction of the sequestration potential. However, the traditional ANN still has the problems of slow convergence speed and low accuracy, etc. Swarm intelligent algorithms have gained the attention of many researchers because of their fast convergence speed and better global search ability, and the lion swarm optimization (LSO) algorithm is a new type of intelligent algorithm that exhibits better global convergence and robustness compared to the particle swarm optimization (PSO) algorithm, grey wolf optimization (GWO) algorithm, etc. [26]. Therefore, in order to overcome the shortcomings of the back propagation (BP) neural network with a slow convergence speed and tendency to fall into a local optimum, this study combines the method of chaotic tent mapping and the difference evolution (DE) mechanism to improve the LSO algorithm and then proposes a CO 2 maximum injection pressure prediction model based on BP neural networks optimized by the improved LSO (ILSO) algorithm. This study was aimed at analyzing the variation characteristics of the formation stress field during CO 2 injection through numerical simulation and constructing a model based on optimized neural networks for predicting the maximum injection pressure so as to improve the calculation efficiency and prediction accuracy. The contributions and innovations of the study are as follows.
(1) For the study area, a hydromechanical coupled model has been constructed using the ABAQUS software based on the ideal elastoplastic constitutive model, the stress field distribution, and Mohr-Coulomb criterion.

Numerical Simulation of Maximum Injection Pressure
ABAQUS is widely used for simulations in engineering fields for its powerful data analysis ability [27,28]. In this study, based on the geological formation and hydrogeological conditions of the study area, the damage changes of the caprock and the maximum injection pressure has been analyzed by ABAQUS. In this section, we present the geomechanical model and numerical simulation.
2.1. Caprock Failure Criterion. Studies have shown that pore pressure has a direct effect on the deformation and damage of porous media [29]. The caprock is an elastoplastic body, which will break and weaken when the load reaches its yield strength, and its failure is controlled by effective stress [30]. When the reservoir is injected with a fluid or a gas, the general reduction in the effective stress might lead to the damage of the caprock and reduce its sealing capacity, thus triggering the risk of gas leakage [31]. When the effective stress is lower than the tensile strength, tensile failure will occur. Prediction of the tensile failure is relatively easy due to the relatively low tensile strength of most sedimentary rocks, which is usually assumed to be zero. In this study, based on the Mohr-Coulomb criterion and the maximum tensile stress theory, we have defined the caprock shear failure threshold, f b , as given by Equation (1), where C and α are the rock cohesion and the internal friction angle, respectively; σ 1 ′ and σ 3 ′ are the maximum and minimum effective principal stresses, respectively. When f b > 1, the caprock belongs to stable state, whereas when f b ≤ 1, shear failure occurs. The smaller the value of f b , the more serious is the damage of the caprock.
In this paper, it is assumed that the reservoir skeleton deformation follows the Terzaghi effective stress principle, as shown in Equation (2), where σ ij ′ is the effective stress, σ ij is the total stress, P is the pore fluid pressure, δ ij is the Kronecker tensor, and β is the Biot coefficient.
The stress balance equation is shown in Equation (3), where F i is the volume component.
The fluid flow in the reservoir is governed by Darcy's law. Taking the influence of gravity and capillary force in the flow process into consideration, the gas-water flow equation can be deduced from the law of mass conservation.
The input parameters of the equations that ABAQUS solves include density, elastic modulus, Poisson's ratio, per-meability, and porosity. The output parameters include pore pressure, stress, and the shear failure threshold of the cap, namely, f b , etc.

Geomechanical Model
2.2.1. Geological Survey and Conceptual Model. Taking a subsurface aquifer in a block as the study subject, fully utilizing the measured data of cores, it is known that the cover layer is thickly bedded green and brown mudstone with a maximum single-layer thickness of 50 m, a stable lateral distribution, and good seal. The reservoir is buried at a depth of 1050 m, with a cumulative thickness of 100 m, dominated by fine sandstone, with an average porosity of 19% and an average permeability of 0:1342 × 10 −6 m 2 , which has good physical properties. Considering the sealing safety, it is necessary to carry out the mechanical integrity analysis of the caprock for the injection process. The numerical model adopts the flow-solid coupling seepage model for single-phase fluids in porous media, and the three-dimensional situation of straight-well injection has been simulated by an axisymmetric method centered on the injection well according to the requirements for numerical calculations. The conceptual model is shown in Figure 3(a).
The constructed model extends vertically from the surface to a depth of 2500 m and includes four typical layers, namely, the overburden, the caprock, the reservoir, and the bedrock. No fault through the caprock is found at the top of the overburden, and the injection point is located in the middle and lower part of the reservoir at the symmetry axis. Based on the delineation criteria for finely portraying the vertical thickness and planar spreading scale of the reservoir, as well as the time and efficiency of the simulation calculation, the model has been meshed using an eight-node hexahedral displacement cell and a pore pressure cell. The grid has been divided equally according to the reservoir thickness in the vertical direction. A dense to sparse grid has been used in the horizontal direction, with a small-sized grid near the injection well and a gradually coarsening large-sized grid in the distant area, as shown in Figure 3(b).

Material Properties and Boundary Conditions.
Without considering the heterogeneity, the reservoir caprock parameters were set to the general values of the study area, as shown in Table 1.
The internal friction angle, the maximum horizontal principal stress, and the minimum horizontal principal stress given in Table 1 are mainly used for the calculation of the f b in Equation (1) to judge whether the shear failure occurs in the cap. The maximum horizontal principal stress and minimum horizontal principal stress were derived from the test results of the differential strain in situ stress. The value of the internal friction angle came from the reference [32]. Ignoring the effect of boreholes, the reservoir was injected with gas at constant pressure, the inner boundary of the caprock was nonpermeable, the outer boundary was a constant pressure boundary, and the normal displacement around the model was fixed.
The detailed initial boundary conditions are as follows. The mechanical boundary is the upper boundary without constraint; the lower boundary is normal displacement con-straint; the left boundary is axially symmetrical, and the right boundary is a normal constraint. The permeable boundaries are as follows. The upper boundary is the seepage boundary; the lower boundary is undrained; the left boundary is undrained, and the right boundary is a closed boundary. The initial conditions also include the initial pore ratio, the initial geostress field, and the initial pore pressure given before the model calculation.

Numerical Simulation of the Maximum Injection
Pressure. Based on the ideal elastic model, a numerical simulation study of gas injection in the reservoir was carried out to analyze the damage change of the caprock under different injection pressures. For the sake of clear presentation, only the visualization results of the analysis area (20 m × 100 m for the reservoir and 10 m × 100 m for the caprock) have been presented.

Simulation Results
Based on ABAQUS. The distribution of the initial ground stress was simulated before gas injection to bring the model to equilibrium. Subsequently,

Geofluids
the reservoir was injected with CO 2 at constant pressure, the simulation time was set to 1 year, and 10000 steps were set to calculate and obtain the pore pressure, horizontal stress, and vertical stress distribution before and after injection, as shown in Figures 4-6, respectively. For a clear representation of the temporal and spatial variations before and after the injection, the variation curve of the effective stress at the measuring point at the reservoir caprock interface with injection time is shown in Figure 7.
Comparing the results of the temporal and spatial changes of the stress field and formation pore pressure in the reservoir caprock from the beginning of the simulation to the end of one year of injection, the following conclusions can be drawn.
(1) The injection of CO 2 has a significant effect on the reservoir pore pressure, which directly affects the effective stress of the formation and the mechanical integrity of the caprock. With CO 2 injection, the formation pore pressure gradually increases, the effective stress significantly decreases, and the maximum stress reduction appears near the reservoir injection zone (2) After injection, the magnitude of vertical stress reduction was greater than that of horizontal stress reduction. One year after the simulated injection, the distribution trend of the pore pressure is consistent with the effective stress of the formation (3) With the increase of injection time, the vertical effective stress and horizontal effective stress at the measuring point gradually decrease, and the change rate of effective stress within one month of injection is large and then gradually decreases and tends to be stable (4) Under the influence of CO 2 injection, the amplitude of the effective stress change at the bottom of the caprock is significantly smaller than that of the reservoir, and the maximum stress change occurs at the intersection of the reservoir and the caprock 2.3.2. Analysis of Reservoir Caprock Damage. Varying the injection pressure at the bottom of the well, the shear failure threshold, f b , defined in Section 2.1, was used to assess the damage to the reservoir caprock during CO 2 injection. The value of f b was represented by SDV5 in the constructed model. Figure 8 shows the SDV5 nephogram of the damage occurred in the reservoir caprock in the equilibrium state before injection and 1 year after the injection at different bottom hole injection pressures. As can be seen from Figure 8, when the injection pressure at the bottom of the well is 22 6 Geofluids the reservoir caprock junction is 1.001, and the caprock is still in a stable state. As the injection pressure increases to 22.530 MPa, the SDV5 reaches 1, and the caprock starts to undergo shear damage. In summary, with the increase of the injection pressure, the caprock gradually began to be damaged from the initial stable state, which could produce the risk of gas leakage. From the simulation results, it is known that the bottom-of-the-well maximum injection pressure for maintaining the integrity of the caprock is 22.529 MPa.

Uncertainty Analysis of Formation Parameters
When using ABAQUS to simulate the maximum injection pressure, the parameters used in the calculation model are empirical or average values of the study area, but measurement and sampling errors and the complexity of the geological environment can lead to some uncertainty in physical and rock mechanics parameters. For example, permeability is related to porosity, and neither has a fixed value but varies within a certain range depending on the actual geological conditions. However, a fixed value is used in the model calculations, which leads to uncertainty in the simulation. In order to obtain results that are closer to the actual geological conditions, the MCS method was used to perform the uncer-tainty analysis of the relevant formation parameters affecting the maximum injection pressure.

Uncertainty Analysis Methods
3.1.1. MCS. MCS is a numerical calculation method based on probabilistic statistical theory, which can accurately analyze the effects of various uncertainty factors to solve the probability of occurrence of a certain event or the expected value of a certain random variable. It is a relatively accurate method for the probability of failure calculation of multidimensional problems, which has been widely used in risk simulation analysis in many fields such as hydropower machinery [33] and geological engineering [34]. In this study, we use Python and MCS to perform an uncertainty analysis of the factors influencing the maximum injection pressure. The analysis contains four basic steps.
(i) Describe the distribution type of parameters and construct a probabilistic model to make the solution correspond to certain characteristics of the random variables of the model   (4). In this study, the elastic modulus and Poisson's ratio obeyed the triangular distribution as known from the field reports and core data of the study area.
If the random variable obeys a Gaussian distribution with an expectation value of μ and variance σ 2 , denoted N ðμ, σ 2 Þ, it is modeled using a normal distribution with the probability density function given in Equation (5). In this study, it is known that the porosity and permeability conform to the normal distribution.   8 Geofluids parameter was obtained by interpreting the logging data, and its probability density function was determined. Further, the influence of different values in the variation range of each parameter on the calculation results of the maximum injection pressure was analyzed by numerical simulation. The relationship between the maximum injection pressure and the formation parameters was established by polynomial fitting. On this basis, the prediction model for the failure probability of caprock was constructed, as shown in In Equation (6), P is the failure probability of caprocks, gðxÞ is a function related to the maximum injection pressure, and f ðxÞ is the probability density function. Since it is difficult to calculate the failure probability using analytical methods, the MCS and the Latin hypercube sampling method were used for dividing each dimension of the 3D model uniformly into m intervals that do not overlap. A point was randomly selected from each interval in each dimension to obtain m 3 subsamples, with each sample point denoted as x j i ðj = 1, 2, ⋯, kÞ, where k is the number of sensitive parameters, and each sensitive parameter in the sample point has its own distribution law. Different sensitive parameters were used as model uncertainties to conduct a large Then, take the mean value as the failure probability of the maximum injection pressure during the initial rupture under different parameters; thus, the distribution model could be constructed. In this paper, the simulation times is 10000 and there are three sensitive parameters; namely, the value of m is set to 100 and k is 3 in our experiment.

Uncertainty Analysis of Model Parameters.
Based on the above model, the uncertainty analysis of the influencing factors of maximum injection pressure was carried out by performing the following steps.
(i) The uncertainty of the physical parameters and rock mechanics parameters of the reservoir-caprock were analyzed, and the sensitive parameters affecting the maximum injection pressure were determined. In this study, we focused on three parameters, namely, permeability, elastic modulus, and Poisson's ratio (ii) The influence law of the uncertainty of different parameters on the maximum injection pressure was analyzed, and polynomial fitting was done using MATLAB. In this study, the fit was characterized by the coefficient of determination (R 2 ) and the RMSE. The closer the R 2 is to 1, the better is the fit; the closer the RMSE is to 0, the better is the fit (iii) The probability density functions of different parameters and the number of simulations were determined, and the MCS-based parameter uncertainty analysis was implemented using Python. In this study, we know that the permeability conforms to a Gaussian distribution, whereas the elastic modulus and the Poisson's ratio obey the triangular distribution (iv) The probability distribution, the cumulative probability graph, and the statistical results were plotted, and the probability distribution interval of the maximum injection pressure based on the influence of different parameters was obtained

Uncertainty Analysis of the Maximum Injection Pressure.
The maximum injection pressure during CO 2 injection is influenced by several factors. This section presents the uncertainty analysis of the formation parameters using MCS and MATLAB based on the actual geological condition.

Distribution of Formation Parameters.
From the logging data and core analysis data, it is known that the third section of the Quantou group of the study area is a "medium-porosity, medium-permeability" reservoir, and its porosity and permeability distribution is shown in Figure 9, with the elastic modulus of 0.93-1.01 GPa and Poisson's ratio of 0.23-0.36.

Analysis of the Effect of Different Parameters on the Maximum Injection Pressure.
To investigate the effect of the parameter uncertainty on the calculation results, the influence of different values on the maximum injection pressure was analyzed using the MCS method within the range of three parameters, namely, permeability (K), elastic modulus (E), and Poisson's ratio (μ). Further, the relationship curves between the maximum injection pressure and the corresponding parameters were plotted (as shown in Figure 10). The maximum injection pressure versus K, E, and μ was obtained by fitting them using MATLAB as shown in Equations (7)-(9), respectively. The R 2 values of the three relations are 0.995, 0.9862, and 0.9956, respectively, and the RMSEs are 0.1147, 0.0073, and 0.0791, respectively, all of which are well-fitted. These relational equations can be used for predicting the maximum injection pressure at different values of K, E, and μ in order to fully assess the risk of caprock damage during injection.  Figure 9: Distribution of porosity and permeability.

Probability Distribution of the Maximum Injection
Pressure. Based on the corresponding fitted relational equations and the probability density functions of K, E, and μ taken randomly, the calculation results with certain probability distribution were obtained by applying MCS sampling. The interpretation of the logging data showed that the reservoir permeability obeys a normal distribution (136.3193 and 137.8410), the elastic modulus obeys a triangular distribution (0.93, 0.97, and 1.01), and Poisson's ratio also obeys a triangular distribution (0.23, 0.29, and 0.36). After determining the probability density functions of the three parameters, the maximum injection pressure was defined as the prediction parameter, and the probability distribution and the cumulative probability distribution of the statistical maximum injection pressure were calculated by simulation using the direct sampling method and setting the number of sampling times to 10000, as shown in Figure 11.
Since several formation parameters have uncertainties, the maximum injection pressure uncertainty analysis should be based on the distribution of each parameter for a compre-hensive assessment of the caprock failure probability. As can be seen from Figure 11, the maximum injection pressure intervals corresponding to the 5% to 95% confidence level for K, E, and μ are (22. 22.529 MPa, which is lower than the value corresponding to the confidence level of 5% for the three parameters, indicating that the possibility of the caprock failure is relatively small. The uncertainty analysis of the formation parameters can provide a more comprehensive and accurate assessment of the maximum injection pressure and provide a powerful guide for the safe injection management of the CO 2 geological sequestration.

BPNN Based on ILSO Algorithm
The method of fixing other parameters and analyzing only the effect of a single factor on the maximum injection pressure has some limitations because of the correlation between certain formation parameters. In contrast, ANNs do not require a prior determination of mathematical equations for mapping relationships between inputs and outputs and have powerful nonlinear approximation capabilities. Thus,

Geofluids
the BPNN was used to predict the maximum injection pressure in this study. To address the problems of the sensitivity of the BPNN to the initial weight and threshold, its slow convergence speed, and tendency to fall into a local optimum, the LSO algorithm was improved by combining the chaotic tent mapping with the DE mechanism, and the ILSO algorithm was applied to optimize the network weights and threshold, in order to improve the global search ability, convergence speed, and prediction accuracy of the model.  [37]. In the LSO algorithm, the search range becomes smaller in the late stage of the optimization for the lion population individuals, thus making the algorithm prone to early convergence and local optimum, which reduces the diversity and the global search ability. To this end, an improved chaotic tent mapping sequence was introduced to initialize the population for its randomness, ergodicity, and regularity which can be used for improving the diversity and uniform ergodicity of the population distribution, inhibit the algorithm from falling into a local optimum, and enhance its global search capability.
Additionally, the mechanism of lioness position update in the traditional LSO algorithm exists only within the lioness population, and the effective information of the lion king   12 Geofluids and cubs might be wasted to a large extent and affect the convergence speed of the algorithm. Hence, it was improved so that when the lioness updates her position, she learns from the lion king and other cubs, in addition to communicating with the lioness group to enhance the global search capability. Further, a DE mechanism was introduced in the lioness position update to improve the accuracy while enhancing the local search capability of the algorithm.

Population Initialization Based on Chaotic Tent
Mapping. A tent map involves segmented linear mapping, which is also a two-dimensional chaotic mapping with a uniform distribution function and good correlation in its parameter range. It is expressed by Equation (10). A tent map has uniformity and randomness, which can make the algorithm easily escape from the local optimal solution so that the diversity of the population can be maintained and the global search ability can be improved [38]. However, there are small cycles and unstable cycle points in its sequence. Thus, a random variable was introduced into the original chaotic tent mapping, and the improved chaotic tent mapping expression is given by Equation (11).
In Equation (11), N is the number of particles within the sequence and rand ð0, 1Þ is a random number ranging between [0,1]. In this study, we use the improved chaos tent mapping to initialize the population. By intro-ducing the random variables, not only does it maintain the randomness, ergodicity, and regularity of the chaotic tent mapping but can also efficiently avoid the iterations from falling into small and unstable periodic points, thus solving the problem of falling into local minima during the training process and consequently strengthening the global search capability of the algorithm.

Lioness Position Update Based on DE. DE algorithm is
an efficient global optimization algorithm with a simple structure, easy implementation, fast convergence, and robustness [39]. The original DE algorithm includes three operations of mutation, hybridization, and selection [40]. This study mainly uses the mutation operation of the DE algorithm to optimize the position update strategy of the lioness. Assuming that x k i denotes the i th particle in the k th iteration, the mechanism of the mutation operation is expressed by Equations (12) and (13), where x k g , x k h , x k j represent the different individuals in the k th generation, F is the scaling factor, and F0 is the variance operator. On constant updating, F gradually converges to F0, which can avoid the loss of the optimal solution of the parent generation and preserve it for the offspring.
In a D-dimensional target search space, a group of N lions is formed, among which only one is a male lion, and the rest are lionesses with a certain proportion of adults. The position of the i th lion is denoted as x i = ð x i1 , x i2 , ⋯, x iD Þ, 1 ≤ i ≤ N. Different types of lions move their positions differently. The lion king moves in a small area at the best food to secure his privilege, and his position is updated according to Equation (14). The lioness requires the collaboration of other lionesses and her position is updated according to Equation (15). The process

F5
Booth 32,32] (0,…,0) 13 Geofluids of the lioness position update communicates only with the lioness group and ignores the valid information of the lion king and cubs. To enhance the global search capability of the algorithm, the lioness position update strategy was improved by using the DE mechanism, as shown in Equation (16).
x k+1 14 Geofluids In Equations (14)-(16), x k+1 i denotes the position of the current lion i in the ðk + 1Þ th iteration; g k is the optimal position in the k th generation; p k i is the historical optimum in the k th iteration of the current lioness i; p k c is the historical optimum in the k th iteration of a lioness c randomly selected in the population of lionesses other than the lioness i; γ is a random number generated based on the normal distribution Nð0, 1Þ; F is a scaling factor with the same meaning as in Equation (13); the probability factor, q, is a uniform random number generated based on the uniform distribution U½0, 1; and α f is the lioness position update perturbation factor, which is defined by In Equation (17), T is the maximum number of iterations; t is the current number of iterations; step = 0:1 × ðhigh − lowÞ, where high and low are the maximum and minimum values, respectively, on each dimension. The perturbation factor, α f , makes the move step of the lioness change continuously with the number of iterations, first surveying with a larger step, followed by a change in the step size from large to small, and finally searching finely with a small step that converges to zero.

Performance Test of ILSO Algorithm.
To verify the performance of the ILSO algorithm, eight benchmark functions were selected for the experiments, as shown in Table 2, where x is the independent variable of the test function, D is the independent variable dimension, Range represents the range of variables, and X min is the global optimal position. The search history of the ILSO algorithm on these 8 functions is shown in Figure 12, which demonstrates the global best position of each function; the historical positions of the lion king, lioness, and lion cub in each iteration; and the convergence curve of the optimal fitness.
Additionally, the performance of the ILSO algorithm was compared with PSO, GWO, and LSO algorithms. In these experiments, the same parameter settings were used for all algorithms; the population size was 50, the maximum number of iterations was 6000, and 50 independent runs on each test function, and the mean value (Mean), standard deviation (Std), and best value (Best) were used as the evaluation metrics. Experimental results are presented in Table 3, where C represents the characteristic of functions, including unimodal separable (US), unimodal nonseparable (UN), multimodal separable (MS), and multimodal nonseparable (MN), and F min is the global optimal value. These eight functions are all representative, among which F5 to F8 are multimodal functions, F7 has several local extremes, and F8 is complex due to its uneven search surface, which makes the search process more complicated.
Experimental results in Table 3 and Figure 12 show that the ILSO algorithm can make the search process more refined and obtain accurate results by introducing the chaotic tent mapping strategy and the DE mechanism in the face of the more complicated search process and the difficulty of finding the optimal value. In addition, when the change of the search value falls into a long-time stagnation, the search value can still continue to fall, which indicates that the improved algorithm has the ability to jump out of the local extremes. Compared to the other three algorithms, the ILSO algorithm proposed in this study has better accuracy and stability.

ILSO_BP Model.
To address the problem of the random initialization of weights and thresholds in BPNN which leads to poor generalization ability and its tendency of falling into local optimal value, the ILSO algorithm was used for optimizing the initial weights and thresholds in order to improve the convergence speed and accuracy.

BPNN.
The BPNN is a multilayer feed forward neural network trained by the back propagation of errors. It uses gradient descent and search techniques to continuously adjust the weights and threshold, which minimizes the MSE between the actual output and the desired output. The method has a strong nonlinear mapping ability, high self-learning, self-adaptivity, and generalization ability. It has been widely used in the field of geological engineering for lithology identification [41], logging interpretation [42], and reservoir prediction [43]. The traditional BPNN is a local search optimization method, and the network weights are gradually adjusted by following the direction of local improvement, which easily leads the algorithm to fall into local extremes, thus leading to the failure of network training. In addition, because the search process of the BPNN uses the gradient descent method, when the objective function to be optimized is very complex, the algorithm will be   16 Geofluids inefficient and the convergence speed is often slow. Therefore, many researchers optimize BPNN by swarm intelligent algorithms in order to improve their convergence speed and global search capability and avoid falling into local optima.

Implementation of ILSO_BP Model.
The process of the ILSO_BP model is shown as Figure 13. The main idea of the ILSO_BP model is to introduce chaotic tent mapping and the DE mechanism in the LSO algorithm to improve the population initialization distribution and the lioness position update method, respectively, and apply the ILSO algorithm to the optimization of the initial weights and thresholds of the BPNN to enhance the global search capability, realize a more extensive and uniform search, and achieve fast convergence speed and high accuracy.

Prediction of the Maximum
Injection Pressure 5.1. Model Parameter Setting. Experience shows that a threelayer BPNN can approximate a function of arbitrary complexity with arbitrary accuracy. Thus, in this study, a threelayer BPNN was used for predicting the maximum injection pressure for CO 2 sequestration, and the initial weights and thresholds of the network were generated and trained by the ILSO_BP model. The input parameters are the four causeless variables affecting the maximum injection pressure, namely, porosity, permeability, elastic modulus, and Poisson's ratio, and the output is the maximum injection pressure. According to Kolmogorov's theorem, the number of nodes in the hidden layer was determined to be four, the transfer function of the intermediate layer neurons was Tanh, the initial learning rate was set to 0.01, the Adam algorithm was used for gradient update, and the number of iteration steps was set to 3000 and 6000, respectively. The population was set to 50, the proportion factor, β, of adult lions was 0.2, and the value interval of weights and thresholds was [-5, 5].

Experimental
Results and Analysis. The 140 sets of sample data obtained from the above ABAQUS simulation were selected, and the training and test sets were divided in the ratio of 8 : 2. MSE, RMSE, and MAE were used as evaluation metrics for model performance, as shown in Equations (18)- (20), where y i is the true value, b y i is the predicted value, and N is the number of samples. The proposed ILSO_BP model was compared with BP, LSO_BP, and GWO_BP methods. The iterative curves of the training errors of the four methods are shown in Figure 14. The final training error and testing error of each method are given in Table 4.

MAE
From Figure 14 and Table 4, it can be seen that the ILSO_ BP model exhibits the lowest training errors as well as testing errors, a higher prediction accuracy compared to the other three algorithms, and the fastest convergence rate.
In addition, the errors between the actual values and the predicted values of the 28 data in the test set were counted, as shown in Figure 15. As can be seen, compared to the other three methods, the experimental results of the ILSO_BP model are all satisfactory, with fewer predicted outliers, and the overall error with the actual value is within a small range. The method proposed uses the ILSO algorithm based on chaotic tent mapping and DE mechanism in the optimization of the BPNN, which solves the defects of slow convergence and low training accuracy and achieves a higher calculation speed and accuracy in prediction of the maximum injection pressure.

Conclusion
A coupled hydromechanical model was constructed using ABAQUS to simulate the mechanical changes of the caprock during CO 2 injection, and the maximum injection pressure was predicted based on the ideal elastoplastic model, Mohr-Coulomb, and tensile failure criteria. Considering the influence of the uncertainty of formation parameters such as porosity and permeability on the reliability of the computational model, the uncertainty analysis of the relevant parameters was carried out using the MCS method to achieve a comprehensive assessment of the probability of caprock damage. Based on simulation results, a prediction model was constructed. The LSO algorithm was improved to optimize the structure of the BPNN by combining the chaotic tent mapping with the DE mechanism to address the problems of slow convergence of the BPNN and its tendency to fall into local optimal solutions. Experimental results show that the ILSO_BP model achieves better performance in the comprehensive evaluation of the maximum injection pressure for CO 2 geological storage, and the   18 Geofluids convergence speed and prediction accuracy are higher than those of other three methods. Thus, the method proposed in this work can provide useful guidance for the injection management of geological CO 2 storage projects, improve safety and storage efficiency, and reduce the risk of gas leakage. However, the actual geological situation is much more complicated, and further research is required to evaluate the integrity and sealing of the caprock during CO 2 injection more comprehensively and accurately. In particular, the presence of faults in the sealing zone during injection will be an important research topic in the future. In addition, improving the model calculation speed further will also be the next work.

Data Availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.