Optimization of Operation Parameters for Helical Flow Cleanout with Supercritical CO2 in Horizontal Wells Using Back-Propagation Artificial Neural Network

Sand production and blockage are common during the drilling and production of horizontal oil and gas wells as a result of formation breakdown. The use of high-pressure rotating jets and annular helical flow is an effective way to enhance horizontal wellbore cleanout. In this paper, we propose the idea of using supercritical CO2 (SC-CO2) as washing fluid in water-sensitive formation. SC-CO2 is manifested to be effective in preventing formation damage and enhancing production rate as drilling fluid, which justifies tis potential in wellbore cleanout. In order to investigate the effectiveness of SC-CO2 helical flow cleanout, we perform the numerical study on the annular flow field, which significantly affects sand cleanout efficiency, of SC-CO2 jets in horizontal wellbore. Based on the field data, the geometry model and mathematical models were built. Then a numerical simulation of the annular helical flow field by SC-CO2 jets was accomplished. The influences of several key parameters were investigated, and SC-CO2 jets were compared to conventional water jets. The results show that flow rate, ambient temperature, jet temperature, and nozzle assemblies play the most important roles on wellbore flow field. Once the difference between ambient temperatures and jet temperatures is kept constant, the wellbore velocity distributions will not change. With increasing lateral nozzle size or decreasing rear/forward nozzle size, suspending ability of SC-CO2 flow improves obviously. A back-propagation artificial neural network (BP-ANN) was successfully employed to match the operation parameters and SC-CO2 flow velocities. A comprehensive model was achieved to optimize the operation parameters according to two strategies: cost-saving strategy and local optimal strategy. This paper can help to understand the distinct characteristics of SC-CO2 flow. And it is the first time that the BP-ANN is introduced to analyze the flow field during wellbore cleanout in horizontal wells.


Introduction
There are many intelligence methods that can be applied in data fitting: Genetic Algorithm (GA) that simulates natural selection and heredity process, Simulated Annealing Algorithm (SAA) that simulates the stochastic state of solid particles cooling from high temperature, Artificial Immune Algorithm (AIA) that is based on the biological immunity theory, Artificial Neural Network (ANN) that simulates the signal transition between neurons, Support Vector Machine (SVM) et al. Each method has its strength and weakness. For example, ANN is very efficient at data fitting but over-fit may happen. SVM is a simple but robust tool for data fitting, especially for small sample and non-linear problems. But it becomes less efficient with large data size [24,25].
There are some recent progresses in fluid flow analysis using complex network that need to be addressed. Complex network has enjoyed a noteworthy development in the last decade and it provides a useful framework to investigate complex systems (in this case, complex fluid flows) from different perspectives [26][27][28][29][30]. To build a complex network from a complex system, system components (flow signals) are represented as nodes while the interactions between nodes are regarded as edges. In a real complex network, there are 'community' structures that incorporate certain group of nudes. Nodes are compactly interconnected within a community and the links between communities are sparse. Community detection is of great significance for clarifying the structure of complex networks. A reliable method for constructing directed weighted complex network from time series data was established in [26]. The method introduced was a faithful tool to extract the dynamical information from experimental signals in complex systems. In order to uncover the transitional behavior of slug flow, Gao et al. [27] calculated multivariate pseudo Wigner distribution (MPWD) and the multivariate multiscale sample entropy (MMSE) for different flow conditions. The results indicated combining the MPWD and MMSE enabled to reveal the transient and multiscale flow behavior from slug flow to churn flow. What's more, charactering the dynamic behavior and flow structure of twophase flow, which is a challenging problem of significant industrial importance, can be accomplished using complex network. Complex networks with various properties corresponding to different situations have been successfully built to investigate this issue. In experimental horizontal oil-water two-phase flow, a complex network-based method was proposed to distinguish complex flow patterns [28]. The results suggested that the community detection in two-phase flow complex network enabled to objectively distinguish intricate horizontal oil-water flow patterns, while the conventional method based on adaptive optimal kernel time-frequency representation (AOK TFR) was invalid. Studies on similar problems have been accomplished with methods based on complex network: the transitions and nonlinear dynamic behaviors of gasliquid flow patterns were quantitatively uncovered based on the distinct topological structures of multivariate complex networks derived from different flow conditions [29]. A new approach based on multi-frequency complex network was proposed to uncover the horizontal oil-water flow structures from experimental multivariate measurements [30]. It was revealed that the community structures could robustly represent the structural features of different flow patterns. In [31], the vertical upward oil-water two-phase flow, which was a multiscale, unstable and non-homogenous complex system, was analyzed using a multivariate multiscale complex network. The results suggested that the clustering coefficient entropy from complex network could not only indicate the oil-in-water flow pattern transition but also show the dynamic behavior of vertical two-phase flow. Based on these applications, the ability of complex network on data analysis of two-phase flow has been fully manifested and it is believed complex network is highly potential to be applied in other fluid flow analysis.
To effectively analyze and fit the data, the BP-ANN was employed in this paper. BP-ANN has been applied with great adaptation and generalization to many other applications in petroleum engineering. Yu et al. [32] combined an ANN with the genetic algorithm and simulated annealing to predict oil reserve quantities based on geological data. Wang et al. built a mathematical model of the main dimensions for self-elevating drilling units by means of BP-ANN [33]. Kaydani et al. and Irani et al. applied an improved ANN to fit core permeability with well logging and geological data [34,35]. The prediction of permeability in a homogeneous area can also be achieved by an ANN with other intelligence algorithms [36]. The bottom hole pressure was predicted precisely without flow pattern determination in underbalanced drilling by applying an ANN [37]. Wang et al. made a 9-variable model based on an ANN to select appropriate deepwater floating platforms [38]. Peng et al. predicted the working life of coiled tubing based on the BP algorithm of an ANN [39]. El-Abbasy et al. developed an experience-based neural network model to predict the rate of failure and working conditions for oil and gas pipelines [40].
Previous cases proved that an ANN shows excellent accuracy when disposing of fitting, regression and prediction problems. Thus, it can be a feasible method for analyzing simulation data and predicting flow field features. On the other hand, SVM also seems to be a qualified method and there have been many comparison studies between ANN and SVM. Though it seems that in most cases SVM has better performance than ANN [41][42][43][44][45][46][47], there are some cases vice versa [48]. Besides, they can also have similar accuracy under certain circumstance [49,50]. Thus, in this paper SVM (more specifically, Support Vector Regression, SVR) was also applied to be compared with BP-ANN in terms of efficiency and accuracy.
The objectives of this paper include are as follows: 1) simulating the SC-CO 2 helical flow field during cleanout; 2) investigating the influence of operation parameters on the SC-CO 2 helical flow field; 3) comparing the differences between SC-CO 2 helical flow and water flow; 4) and optimizing the operation parameters employing BP-ANN approach based on the cost-saving strategy (CSS) and local optimal strategy (LOS).  nozzles and rear nozzles. Each group of nozzles can produce high-pressure SC-CO 2 jets during circulation and has specific functions. Impinging forward, the forward jets are designed to break up encountered consolidated deposits into floating particles. The lateral jets are not orthogonal to the casing. Instead, they impact the casing at predefined angles, thus generating a strong helical flow in the annulus. The rear jets can improve the sweeping efficiency by pushing floating solids backward.

Material and Methods
The concentric annulus and well bottom are regular geometries where flow field is relatively stable, so they were meshed by structural grids. In contrast, the flow field is much chaotic near the exit regions of nozzles. Based on this fact, the whole model was divided and meshed in three sections, as displayed in Fig 3. Particularly, to characterize the flow field precisely and accelerate convergence, nine infill frustum grids were added along the initial jet paths. Corresponding boundary types of walls, velocity inlet, pressure outlet and interfaces were defined.
2.1.2 Mathematical models. Governing equations. Continuity equation where v i (i = x, y, z) is the velocity component in the i direction; ρ is fluid density; t is time.  Momentum equation where g i is the component of gravity acceleration in i direction; @x i is the partial differential in i direction; μ e is effective viscosity; R i is distributed resistance in i direction; T i is viscous loss term in i direction. Energy equation where C p is specific heat; T 0 is total (or stagnation) temperature; K is thermal conductivity; W v is the viscous work term; Q v is the volumetric heat source; F is the viscous heat generation term; E k is kinetic energy; Turbulence model. According to the helical characteristics of the flow field, the RNG k−ε model was chosen to describe the flow field for better performance [6,51], which has a similar form to the standard k−ε model: where G k is the generation of turbulence kinetic energy due to the mean velocity gradients; G b is the generation of turbulence kinetic energy due to buoyancy; Y M is the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate; α k and α ε are the inverse effective Prandtl numbers for k and ε; S k and S ε are source terms. Sate equations. When flowing in annulus, the properties of SC-CO 2 are controlled by temperature and pressure. Thus, the state equations of carbon dioxide should also be considered. In the numerical simulation, the state equation developed by Span and Wagner [52], which is recommended by NIST, was used to calculate the density and isobaric heat capacity of SC-CO 2 . Based on Helmholtz energy, Span-Wagner state equation is a highly accurate reference equation for the thermodynamic properties of pure carbon dioxide. It can be applied with temperature ranging from the triple-point temperature to 1000 K and a maximum pressure of 2200 MPa.
The Helmholtz energy (A) is the function of two independent variables: density (ρ) and temperature (T), namely, A = A(ρ, T). Span-Wagner state equation is a function of dimensionless Helmholtz energy (α), which is given as follows [52]: where δ is dimensionless density (δ = ρ/ρ c ); τ is inverse dimensionless temperature (τ = T c /T); ρ is the density of CO 2 , kg/m 3 ; T is temperature, K; R is the gas constant, kJ/ (kgÁK); T c is the critical temperature of CO 2 , K; ρ c is the critical density of CO 2 , kg/m 3 . The dimensionless Helmholtz energy is composed of two parts: aðd; tÞ ¼ a o ðd; tÞ þ a r ðd; tÞ ð 7Þ where α o is the ideal part of dimensionless Helmholtz energy; α r is the residual part. The density and isobaric heat capacity (C p ) of CO 2 are calculated by equations as follows: where p is pressure, kPa.
In addition, the viscosity and thermal conductivity of CO 2 are calculated using the model presented by Fenghour et al, recommended by NIST [53,54]. In this model, viscosity (η) is divided into two different terms as in Eq 10: where η is the viscosity of CO 2 , PaÁs; η 0 is the viscosity of dilute gas, PaÁs; η R is the residual part of viscosity, PaÁs. The thermal conductivity (λ) is calculated in a similar manner: where λ is the thermal conductivity of CO 2 , W/ (mÁK); λ 0 is the thermal conductivity of dilute gas, W/ (mÁK); λ R is the residual part of thermal conductivity, W/(mÁK); λ c is the critical enhancement of thermal conductivity, W/ (mÁK).

Simulation parameters.
Five main factors influencing the helical flow field were investigated: nozzle assembly, flow rate, ambient pressure, ambient temperature and jet temperature, as shown in Tables 1 and 2. Ten possible assemblies were applied in the simulation with flow rates ranging from 10 L/s to 20 L/s. The ambient pressure/temperature indicates the true vertical depth (TVD) of the well. In this paper, the ambient pressure changes from 10 MPa to 50 MPa, roughly corresponding to a TVD variation from 1000 m to 5000 m. The jet temperature changes from 353 K to 433 K, and the ambient temperature changes from 333 K to 433 K. Both temperature and pressure are determinants of the physical properties of SC-CO 2 . The leakage of CO 2 into the formation is ignored. The total area of nozzles decides the exit velocity of fluid, larger area means smaller exit velocity. The enlargement of a specific type of nozzle will increase its share of total flow rate but decrease the exit velocity for all nozzles since total nozzle area is increased. 2.1.4 Simulation setup. The steady flow situation was defined in all simulations. Governing equations were discretized using second order upwind scheme. The SIMPLE coupling method was used to solve the equations. To improve computation accuracy, the criterion of convergence was set as that all residuals fell blew 1×10 −5 instead of default value 1×10 −3 . The simulations were accomplished using commercial CFD software ANSYS Fluent 14.5. A diagram depicting the whole simulation process is provided as Fig 4. It should be pointed out that when computation is not converged refining mesh and/or changing numerical method will be helpful to improve convergence.

BP-ANN Model
A sophisticated artificial intelligence algorithm, a BP-ANN, is basically a large class of parallel processing structures that are able to simulate vague and complicated connections between inputs and target data through the application of many nonlinear processing units called 'neurons' [55,56]. The connection between inputs and target data can be 'learned' by the neural network after adequate training [57,58]. A three-layer feed-forward neural network with a backpropagation algorithm can map any nonlinear relationship [59]. A potential problem that should be noted when using this powerful nonlinear regression method is over-fit. It is possible that there are some unqualified samples in the data. Over-fit is the match between the erroneous inputs and target data that introduces wrong relationships into the model. It will misunderstand 'noise data' as correct ones, thus impairing the certainty and precision of the connections between inputs and targets.
In this paper, the tangential velocity (V t ) and the annular velocity (V a ) of a SC-CO 2 helical flow field are set to be functions of flow rate, nozzle assembly, ambient pressure, ambient temperature, jet temperature, cleaning distance and radial position (dimensionless radius). A three-layer BP-ANN ( Fig 5) is developed to match the operation parameters with V t and V a , in which the transfer functions in the hidden layer are sigmoid, whereas those in the output layer are linear functions. The input layer consists of 7 neurons that represent the 7 operation parameters affecting the helical flow field, whereas the output layer has 2 neurons that  The simulations cover all the ten potential nozzle assemblies and a large range of parameters in practice.
represent V t and V a . The number of neurons in the hidden layer determines how well a problem can be learned. If there are too few neurons, the network will be more generalized than necessary and fail to learn the specific patterns very well. Otherwise, if there are too many neurons, the network will have the tendency to memorize the specific problem and not be generalized enough for application [51]. In other words, the number of neurons in the hidden layer relies on the nature of the problem and available data. Neither too many nor too few neurons in the hidden layer should be used. The neuron number in hidden layer was set as 10 as the result of considerations of accuracy, training time and the risk of over-fit. To determine the optimal number of neurons in hidden layer, firstly we found a rough range of reasonable number of neurons based on references where the problems are nonlinear and have similar data structures (5~8 inputs, 1~2 outputs compared with 7 inputs and 2 outputs in this paper) [32,[34][35][36]38,40]. Thus, the test range of neuron number was decided as 3 to 14. Then we tried each number and obtained corresponding training time and accuracy. When neuron number was small (3 to 6), convergence was fast but the accuracy was not satisfactory (R 2 were below 0.96). With the increase of neuron number (7 to 11), BP-ANN became more accurate (R 2 were above 0.98) but need more training time. An interesting phenomenon is that after neuron number was raised to 12, the times for iteration and single convergence dramatically increased. In this case over-fit was assumed to occur. So we narrowed down the range of neuron number as 7 to 11. In these range 10-neurons was chosen because it had the second highest accuracy and 11-neurons case with the highest accuracy was right at the edge of over-fit. There may be some deviation of the optimal neuron number. However it is revealed in the figure that the accuracy of BP-ANN model will only have tiny changes if the number of neurons varies a little. Therefore, after several trials, we found that the best number of neurons in the hidden layer is found to be 10. In this way, a BP-ANN with a 7-10-2 topology is established (Fig 5).
During the training process of BP-ANN, the error is subsequently backward propagated through the network to modify the model parameters, such as the weights of different layers and neuron thresholds, by means of the gradient descent method. The purpose is to minimize the sum of the mean squared error (MSE) between network outputs and the target data. In this case, the target data are the real values of V t and V a from CFD simulations. The MSE is given by: where U is the sum of the MSE between the target data and the network outputs, N is the number of input samples, M is the number of output neurons, T j (i) is the ith component target value corresponding to the jth input, and Y j (i) is the network output that approximates the target value T j (i).
In this model, to avoid over-fit, all data are divided into three sets: training dataset, validation dataset and test dataset. During the learning process, the model is only fed the training data. The testing dataset is used to test the network during training and also to correct it continuously by adjusting the parameters (weights, thresholds, etc.). The validation dataset, which is not presented to the network during training, is used to validate the model. The model will be retrained if the validation performance is poor. To perform data classification, each dataset is randomly selected from the entire database according to the classification percentages. In this paper, the percentages of training, validation, and test datasets are fixed to be 70%, 15%, and 15%, respectively.
The model framework is shown in Fig 6. Training iteration will end if a certain number of training epochs are conducted without a further decrease of absolute error between outputs and target data. When this criterion is met, training is stopped. If the validation performance is bad, the transfer functions or the number of neurons in the hidden layer should be changed.

Characteristics of the SC-CO 2 helical flow field
In this part, the outcomes of the CFD simulation are analyzed mainly in terms of two major parameters of the SC-CO 2 helical flow field-i.e., both tangential velocity (V t ) and axial velocity (V a ). In horizontal well cleanout, V t and V a play a significant role in suspending and transporting sands. V t implies the shear force applied to solid particles in the helical flow field, which helps suspend the particle. In general, a large V t will effectively prevent sands from settling down. Similarly, V a implies the normal stress (pushing force) applied to solid particles. The higher V a is, the faster sands move back into the flow field.
3.1.1 General features of the SC-CO2 helical flow field. During the cleanout operation, high-speed SC-CO 2 jets constantly impact the wellbore, forming the annular helical flow field. Pathlines near the cleaning bit and of separate jets are displayed in Fig 6. The lateral jets impinge on the casing slantingly and generate an asymmetry overflow near the casing, which is heterogeneous and evolves to a swirling flow. Then, the initial helical flow mingles with the forward and rear jets in annulus. Based on flow pattern, the whole flow field can be divided into two zones: transition zone and laminar zone. As already mentioned, the flow field is quite chaotic and unsteady near the exit regions of the nozzles, with high turbulence (Fig 7(a)) caused by the impact of rear and forward jets on the casing and the mixing of three jets. We define the zone extending from the cleanout bit to the place where the velocity becomes annularly symmetrical as the transition zone or turbulent zone. After a certain distance from the bit, the flow field stabilizes and finally becomes uniformly helical (Fig 8(b)), which is referred as the laminar zone. Despite the different simulation parameters, we generally find that the length of the transition zone is approximately 0.6 m (Fig 8).
Technically, the SC-CO 2 helical flow field is more turbulent than the water helical flow field in transition zone. With the same pumping rate, temperature and ambient pressure, the velocity and density of SC-CO 2 and water are comparable (for example, the average density of SC-CO 2 under the condition of Fig 9 is 665 kg/m 3 , while water is about 1000 kg/m 3 ). On the other hand, the viscosity of water is over one magnitude larger than that of SC-CO 2 (the viscosity of water is 0.2838 mPaÁs while that of SC-CO 2 is only 0.0278 mPaÁs under the condition of Optimization of SC-CO2 Cleanout Using BP- ANN  Fig 9). So from Re = ρvd/μ we can conclude that the SC-CO 2 flow has higher Reynolds number (turbulence level) than water flow (In the same wellbore model, the characteristic length d is the same for water and SC-CO2 flow). Higher turbulence in transition zone is welcomed because it means there is more energy to stir up consolidated sands.
In laminar zone, the low viscosity of SC-CO 2 leads to an obvious decrease of V t (Fig 9) and a more uniform distribution of V a (Fig 10) in cross sections. The attenuation of V t is a disadvantage in terms of suspending sands, but the uniform distribution of V a means SC-CO 2 is better at sweeping sands. The uniform V a in cross sections can induce a uniform backward movement of sands, which is highly preferred in practice.  Optimization of SC-CO2 Cleanout Using BP-ANN Both V t and V a are close to zero near the casing and cubing, so an effective flow zone is defined as the radial ring-like region whose dimensionless radius is between 0.1 and 0.95. Then, V t,avg and V a,avg can be calculated as the circular-area averaged tangential/axial velocities (Fig 11) that represent V t and V a in different cross sections. The cleaning distance is defined as the axial distance from a certain cross section to the cleanout bit. With regard to V a,avg , there is little difference between the SC-CO 2 helical flow field and the water helical flow field, whereas the V t,avg of SC-CO 2 flow is significantly smaller than that of water flow. Again, this is because SC-CO 2 has a much smaller viscosity than water, while the densities of SC-CO 2 and water are similar.
3.1.2 Effect of different nozzle assemblies on V a and V t . The nozzle assembly, meaning the orifice diameters of each nozzle group, affects the hydraulic power of each high pressure SC-CO 2 jet, which consequently influences the strength of annular SC-CO 2 helical flow.
Lateral nozzle size. Fig 12 exhibits the V t and V t,avg distributions with the same simulation parameters except lateral nozzle size. There is no monotonic relation between V t /V t,avg and the lateral nozzle size. The No.4 nozzle assembly, with 4 mm lateral nozzles, has almost the same V t,avg distribution as the No.8 nozzle assembly, whose lateral nozzles are 6 mm. Lateral jets are the main source power of rotational flow and thus determine V t . With larger lateral nozzles, lateral jets will have larger flow rate, which means that the flow field will be more rotational,  Optimization of SC-CO2 Cleanout Using BP-ANN and larger V t will be achieved. Nevertheless, at the same time the nozzle equivalent diameter will be larger, leading to a reduction of lateral jets exit velocity. In this way, V t is the outcome of the balance between larger lateral jets flow rate and lower exit velocity. In other words, there should be an optimal size of the lateral nozzle.
Rear nozzle size. In contrast with the lateral nozzle, the rear nozzle size has an evident monotonic relation with V t . Fig 13 presents the radial distribution of V t and axial distribution  of V t,avg varying with rear nozzle size. As the rear nozzle size increases, lateral jets have less flow rate, and V t decreases significantly, reflecting the decreased helical flow strength. This decrease of V t is most obvious when the rear nozzle size increases from 5 mm to 6 mm. Additionally, as the cleaning distance increases, the influence of the rear nozzle size on the helical flow strength weakens. The largest difference in V t,avg occurs between the No.1 and No.10 assemblies at a cleaning distance of 0.6 m. V t,avg of the No.1 assembly is 1.81 m/s, whereas that of the No.10 assembly is 1.05 m/s. This can lead to a distinct difference in cleanout efficiency. Minimizing the rear nozzle size seems a possible way to improve the suspending ability of the helical flow field. However, if this is the case, the sweeping efficiency will be damped as the distributed flow rate for rear nozzles decreases. Actually, the rear jets can weaken V t but enhance V a in the annular helical flow.
Forward nozzle size. Compared to the rear nozzle, the forward nozzle size has less influence on V t distributions (Fig 14). As forward nozzle size increases, V t,avg has a much lower tendency to decrease. Likely, minimizing forward nozzle size for the purpose of larger V t is not the best way because the forward jets require enough flow rate and energy to break up solid sediments in the front of the cleaning bit. In conclusion, reducing the rear/forward nozzle size or increasing the lateral nozzle size can enhance the strength of annular helical flow by increasing V t . However, the jets for penetration and sweeping will be reduced. Thus, there exists an optimal nozzle assembly that can provide the best combination of V t and V a . This issue will be addressed later in the BP-ANN optimization model.
3.1.3 Effect of flow rate on V a and V t . Naturally, a higher flow rate will induce larger V t and V a . We find that during the simulation, both V t and V a are almost directly proportional to the flow rate, as are V t,avg and V a,avg (Fig 15). In this sense it seems that the flow rate should be as high as possible. However, a high flow rate requires greater pump and increases energy costs. To strike an economical balance between velocities and flow rate cost, BP-ANN optimization model can be applied; this will be discussed in section 3.2.
3.1.4 Effect of ambient pressure on V a and V t . Ambient pressure is found to be the least influential factor. In fact, V t,avg and V a,avg hardly change with different ambient pressures (Fig  16). Considering the relationship between ambient pressure and TVD, this means that the SC-CO 2 helical flow cleanout can be performed in horizontal wells whose vertical depth is no more than 5000 m, assuming that the pressure gradient is 1 MPa/100 m.
3.1.5 Effect of jet and ambient temperature on V a and V t . Temperature is a minor factor in water helical flow cleanout because the properties of water hardly change with temperature. In contrast, the physical properties of SC-CO 2 are sensitive to temperature. Ambient temperature is defined as the temperature of the formation, casing and tubing, and jet temperature is the temperature of jets when they squirt from the nozzles.
Ambient temperature has an approximately linear relationship with V t,avg (Fig 17) and V a,avg at a certain cleaning distance. V t,avg increases with higher ambient temperature. At cleaning distances of 0.6 m and 2.0 m, the percentages of V t,avg with ambient temperature increasing from 333 K to 433 K, increase by 16.8% and 25.7%. Furthermore, the percentages of V a,avg increase by 26.9% and 32.0%. The increased ambient temperature has a greater influence on V a,avg than V t,avg . Moreover, at larger cleaning distances, the influence of ambient temperature is enhanced.
Analogously, V t,avg and V a,avg also have proportional relationships with jet temperature (Fig  18). However, the difference is that V a,avg and V t,avg decrease with increasing jet temperature. At a cleaning distance of 0.6 m, with the jet temperature increasing from 353 K to 433 K, the percentages of V t,avg and V a,avg decrease by 13.6% and 20.1%. At a cleaning distance 2.0 m, the percentages of V t,avg and V a,avg decrease by 21.1% and 24.1%. Again, V a,avg is more sensitive to jet temperature than V t,avg , and the velocities are influenced more by jet temperature at larger cleaning distances.
Actually, ambient temperature and jet temperature share more common properties. The change of ambient/jet temperature is, substantially, no more than the change in the difference between these two temperatures. We picked cases that had the same difference between the ambient and jet temperatures and then compared the distributions of V t and V a . In the two simulation cases where the ambient/jet temperatures are 373 K/413 K and 333 K/373 K-i.e., the difference between the ambient and jet temperatures is -40K-the distributions of V a , V t , V t,avg and V a,avg are nearly identical (Fig 19). The same phenomenon occurs in all other cases with the same temperature difference. We can conclude that V t and V a depend on the difference between the ambient and jet temperatures, not on these two variables independently.
In SC-CO 2 drilling, liquid CO 2 is pumped into the drilling string. It is heated and pressurized as CO 2 goes deep underground. When the critical temperature and pressure are reached, CO 2 achieves a supercritical state. The relationship between the surface temperature of liquid CO 2 and the temperature of SC-CO 2 jets can be clarified by wellbore heat transfer models [60,61]. Thus, it is feasible to modify the surface temperature based on the design of the jet temperature.
In this section, the features of SC-CO 2 helical flow are investigated and compared with water helical flow. SC-CO 2 helical flow is feasible in specific cleanout cases such as water-sensitive and low-permeability formations. In terms of nozzle assembly, there is an optimal size for the lateral nozzle to produce the best V t . Though a smaller rear/forward nozzle size is beneficial Optimization of SC-CO2 Cleanout Using BP-ANN to improve V t , it may not guarantee the ultimate performance of cleanout. Better flow rate is able to boost V t and V a but requires extra pump and energy costs. Ambient pressure has little influence on SC-CO 2 helical flow. V t,avg and V a,avg improve with higher ambient temperature or lower jet temperature. V t and V a depend on the difference between the ambient and jet temperatures, not on these two variables independently.

Performance of the BP-ANN model
3.2.1 BP-ANN data fitting. In total, 4563 samples covering 28 cases were inputted to the BP-ANN. The training, test and validation processes were performed successfully with sufficient accuracy. There were 212 iterations in all, with a decent rate of convergence, and the best validation performance of MSE was 0.0059 (Fig 20). The error histogram (Fig 21) shows that most absolute errors are within a small range near zero, which is a sign of good performance for the BP-ANN. Moreover, the R 2 values of the training, test, validation and combination thereof are quite close to 1 (Fig 22), indicating perfect regression between the outputs and targets. It should be noted that in regression plots, there are several points that are far away from the diagonal (regression line). These points are abovementioned 'noise' points resulting from either improper sampling position or perturbance of the helical flow field. In this model, these 'noise data' are not included, thus avoiding the over-fit problem.
3.2.2 Comparison between SVM and ANN. SVR (Support Vector Regression) based on SVM method was applied to fit the operation parameters with V t and V a . RBF (Radial Basis Function) was chosen as the kernel function. Since one E-SVR can only have one output, two Optimization of SC-CO2 Cleanout Using BP-ANN independent E-SVRs were trained to match V t and V a , respectively. Similarly as training BP-ANN, total data was divided into two parts: 70% was used for training while 30% was used for validation.
BP-ANN and E-SVR were compared in terms of time consumption and accuracy. The training time for two E-SVRs combined was 8.5s while that for BP-ANN was 23s. Although it seemed that E-SVRs required less time to train, it was finding the best values of g (gamma) for RBF kernel function and c for loss function that costed massive time. Two methods were used to find the best g and c based on cross validation: Grid Method (GM) which consumed nearly an hour; PSO (Particle Swarm Optimization) method which consumed more than 16 hours in total. The sets of 'best' g and c were found to be 256 and 16 for both V t and V a by GM. By PSO, best gs were found to be 100 for V t and V a while cs were 31.8885 for V t and 36.0296 for V a . E-SVR trained with PSO parameters had better accuracy (total R 2 was 0.9530 for PSO-trained E-SVR while that was 0.9439 for GM-trained E-SVR). So E-SVR with parameters by PSO seemed to be the best SVR model for this problem, which had R 2 of 0.9480 and 0.9580 for training and validation. The best MSE of E-SVR was 0.1039 while that of BP-ANN was 0.0059. The R 2 for BP-ANN training, validation and test were 0.9871, 0.9886 and 0.9874 with total R 2 equal to 0.9874. Thus it can be concluded BP-ANN had better performance for this problem than SVM method, both in terms of time (23s vs 16 hours plus) and accuracy (0.9530 vs 0.9874 of R 2 ). And that's why we used BP-ANN instead of SVM.

Optimization strategies.
Once the best match of BP-ANN model is achieved, it can produce predicted V t and V a based on input operation parameters. Combined with corresponding requirements and strategies, optimization of operation parameters can be achieved. Cost-Saving Strategy (CSS). Flow rate is one of the main factors that determine the cost of cleanout. Maintaining a large flow rate requires heavy-duty pumps and excessive energy, which are often expensive. To fulfill the cleanout task, it is only necessary to make sure that the intensity of flow field within a critical cleaning distance is big enough to carry sands. In this way, the essence of CSS is to find the minimal flow rate (pumping rate) that keeps V t,min and V a,min (the minimum tangential and axial velocity) within the effective flow zone of critical cleaning distance above critical values (V tc and V ac , below which sands cannot be carried sufficiently), so that sands can be suspended and carried by the helical flow field. The critical values for cleaning distance, V t and V a are functions of wellbore geometry, flushing fluid property and sand size, which is beyond the scope of present study. V ac sets the upper boundary of the cleanout time and should be decided according to the design requirement. After determining V tc and V ac , we can increase the input flow rate from 10 L/s gradually and get the output increasing V t,min and V a,min from BP-ANN model with other input parameters fixed. The minimal flow rate required is obtained when V t,min and V a,min reach V tc and V ac . The process of parameters optimization using CSS is illustrated in Fig 23(a).
Local Optimal Strategy (LOS). Practically, it is conceivable that only limited types of tools and equipment are available in well sites. In this situation, although the optimal operation parameter set that can have the largest V t and V a cannot be obtained, it is still beneficial if a local optimal set of parameters can be decided. One example is to determine the optimal nozzle assembly when there are only limited types of pumps and cooling equipment available, that is, Optimization of SC-CO2 Cleanout Using BP-ANN when the flow rate and jet temperature have limited options. On the one hand, there is the requirement that V tc and V a,c should be satisfied in the first place; on the other hand, a stronger helical flow field will definitely result in better cleanout performance. In other words, the aims of LOS are: 1) find nozzle assemblies that can meet the requirement of V tc and V ac ; 2) find the best nozzle assembly by which V t,avg and V a,avg in the effective zone can be maximized. The process of parameters optimization using CSS is illustrated in Fig 23(b).
3.2.4 Optimization cases. Case 1. CSS. Assume that there is a horizontal well that requires SC-CO 2 helical flow cleanout. The TVD of the well is 3500 m; the geothermal gradient is 2°C/ 100 m; the surface temperature is 20°C; the pressure gradient is 1 MPa/100 m; the No.3 nozzle assembly is used; the jet temperature is 373 K; the critical cleaning distance is 2.5 m; and V tc and V ac are 0.3 m/s and 0.9 m/s, respectively. Recall the BP-ANN model; it is easy to find the minimum flow rate that meets the requirement of the critical velocities is 15.4 L/s. The corresponding minimum pump rate can thus be selected after other hydraulic losses are considered. In the end, the cleanout can be performed economically with the lowest pump and energy cost.
The V t,avg , V a,avg , V t,max and V a,max distributions along the wellbore are shown in Fig 24. V a,max decreases from the cleaning distance of 0.6 m to 2.0 m, which agrees with the simulation results. However, V a,max has a slight increase when the cleaning distance is larger than 2.0 m. This may be because no sample data have a cleaning distance larger than 2.0 m. The velocities beyond 2.0 m are derived only through interpolation. Nevertheless, apart from V a,max , the other three curves have reasonable good decreasing shapes. The outcome may well be inaccurate if the input has some variables that are far from the simulation ranges. So it is suggested to use the model with inputs whose variables are within the simulation ranges or, less restrictively, not too far from the simulation ranges.
In addition, the other operation parameters such as nozzle assembly and jet temperature, if not fixed, can be optimized as well. The outcome will be the optimal set of operation parameters based on CSS. In the same case, the best nozzle assembly is No.1, and the best jet temperature is 353 K, with a minimal flow rate equal to 14.3 L/s. Compared to the previous result, the optimization of the nozzle assembly and jet temperature successfully reduces the minimum requirement of the flow rate from 15.4 L/s to 14.3 L/s, which can further reduce the cost of the pump and energy. It is proven that the BP-ANN optimization model is valid for optimizing single and multiple operation parameters based on CSS.
Case 2. LOS. Suppose that for a horizontal well that requires a SC-CO 2 helical cleanout, the TVD is 2300 m; the geothermal gradient is 1.6°C/100 m; the surface temperature is 27°C; the pressure gradient is 1.1 MPa/100 m; the jet temperature is 360 K; the cleaning distance is 1.7 m; a 12 L/s pump is the only pump available; all 10 nozzle assemblies are available; V tc is 0.3 m/s; and V ac is 0.8 m/s. Hydraulic losses are ignored for simplification. By recalling the BP-ANN model and programming, the V t,min , V a,min , V t,avg and V a,avg of each nozzle assembly at a values that are larger than V tc , so they are the suitable nozzle assemblies that can satisfy the critical velocity requirement. Furthermore, it can be confirmed that the best nozzle assembly is the No.1 assembly. Optimization of SC-CO2 Cleanout Using BP-ANN Parameter optimizations using the BP-ANN model based on CSS and LOS are successfully performed. The constraint of V tc and V ac plays a key role in determining applicable operation parameters. The two strategies discussed are important but are not the only applications of the BP-ANN model. Another potential application of the BP-ANN model is to find the optimal traveling speed of the tool when performing wiper tripping in cleanout. The application can be further extended to cleanouts in inclined and vertical wells with various flushing fluids if cases that consider deviation angle are provided.

Conclusion
In this paper, CFD simulation is successfully applied to analyze the characteristics of SC-CO 2 helical flow and sand sweeping efficiency in terms of V t and V a . Comparisons are made between SC-CO 2 helical flow and water helical flow. The influences of nozzle assembly, flow rate, ambient pressure, ambient temperature and jet temperature are obtained. The simulation data are then processed by BP-ANN, and the outcome model is able to optimize operation parameters in different situations. In all, the following conclusions can be drawn.
1. Compared to water helical flow cleanout, SC-CO 2 helical flow cleanout, which can stimulate production, is more suitable to be applied in water-sensitive, low-permeability and conventional reservoirs.
2. In a horizontal annulus, the SC-CO 2 helical flow field has more evenly radial distributions of V a and smaller V t due to its low viscosity. The nozzle assembly plays a major role in determining the characteristics of SC-CO 2 helical flow. If other operation parameters are fixed, the suspending ability of SC-CO 2 helical flow obviously improves with increasing lateral nozzle size and decreasing rear and forward nozzle size. For the simulated situations in this paper, the optimal nozzle assembly is identified as the No.1 assembly, whose rear, lateral and forward nozzles are 3 mm, 4 mm and 4 mm, respectively.
3. At certain cleaning distances, V t and V a improve linearly with higher ambient temperature or lower jet temperature. A significant discovery on the influence of temperature is that V t and V a depend on the difference between the ambient and jet temperatures, not on these two variables independently.
4. BP-ANN is successfully applied to match the operation parameters with V t and V a at various axial and radial positions. The outcome model is quite accurate, with R 2 almost equal to 1. SVM is less efficient and accurate than BP-ANN for this problem. It is able to predict V t and V a precisely even if new combinations of operation parameters are provided. However, it is suggested that the model be used with input parameters that are within the simulation ranges, or at least not too far from the simulation ranges, for the purpose of veracity.
5. It is found that the BP-ANN optimization model is valid for optimizing single and multiple operation parameters. Two strategies are presented: the cost-saving strategy and the local optimal strategy. For the two hypothetical cases, optimizations based on the CSS and LOS strategies are performed successfully, which indicates the robustness of the BP-ANN optimization model. There is enough evidence to claim that the BP-ANN optimization model will work very well with field data.
Supporting Information S1 Dataset. Axial and tangential velocities with respect to different ambient pressures.