Multiphase non equilibrium pipe flow behaviors in the solid fluidization exploitation of marine natural gas hydrate reservoir

In recent years, marine hydrate has attracted more and more attentions. As a creative way to efficiently exploit marine hydrate, the method of solid fluidization exploitation is to excavate and crush hydrate, transport hydrate particles to sea surface platform through airtight pipeline, and obtain methane gas finally. When hydrate particles are transported up, as temperature rises and pressure drops, hydrate will decompose at a certain critical position and produce a large amount of gas, bringing complicated gas‐liquid‐solid multiphase non equilibrium pipe flow, and seriously threating well control securities. Thus, we establish mathematical models of wellbore temperature and pressure, hydrate phase equilibrium, hydrate dynamic decomposition in riser pipe flow and wellbore multiphase flow coupling hydrate dynamic decomposition, and propose and verify numerical calculation method. Then numerical model application analysis under different construction parameters is carried out, influence behaviors of multiphase non equilibrium pipe flow are obtained, and field construction measures are put forward. Properly increasing solid throughput can increase production of natural gas. However, problems such as well control risks will intensify. Therefore, delivery capacity, liquid density, and wellhead back pressure should be increased appropriately at the same time. This study provides important technical support for solid fluidization exploitation of marine hydrate.


| INTRODUCTION
The reserves of global natural gas hydrate are huge, which mainly distributed in deep water and land tundra. 1 The reserves of marine gas hydrates are about 100 times as much as that in land. The total conservative estimate of it is 2.83 × 10 15 m 3 , which is 1.56 times as much as the total known natural gas reserves. 2 In the face of such a huge amount of resources, it is of great significance to accelerate the safe and efficient exploitation of natural gas hydrate in the development of world energy.
At present, the samples of marine natural gas hydrate obtained have the characteristics of shallow burial depth and weak cementing between particles. [3][4][5] This type of hydrate is stored in the shallow surface sediments of the deep seabed or is dispersed in the sea floor with particles or veins. The pressure and temperature conditions of these sites can keep hydrate structures stable. However, it does not have a | 761 WEI Et al.
good cap and stable storage space like conventional oil and gas fields, the hydrate itself is the reservoir skeleton instead. It is a kind of metastable and weakly cemented hydrate. [3][4][5] When temperature and pressure in the ocean floor change under the influence of field operation or other ocean engineering, this kind of metastable hydrate will decompose in a large amount, the reservoir will collapse, and the uncontrollable decomposition and dissipation will occur. The occurrence of this situation not only wastes a large amount of resources, but also causes the risk of marine environmental damage, submarine geology, and equipment instability. [3][4][5] Therefore, how to rationally and safely develop the natural gas hydrate in the shallow surface of the seabed is very important for the full utilization of such resources and the development of deep water resources.
The existing exploitation methods of natural gas hydrates are mainly used through reference to the conventional oil and gas reservoir exploitation process. Through the methods of depressurization, heat shock, chemical injection, and replacement, methane gas is released by hydrates decomposition or replacement, and then is produced. Existing examples of exploiting natural gas hydrate reservoir are mainly included as follows: In the Messoyakha gas field of Siberia area, the gas is produced from hydrate decomposition by means of depressurization and chemical injection. [6][7][8] In the Mallik area of Canada, two hydrate samples were tested by means of depressurization and heat injection. 9,10 In the north slope of Alaska, a variety of methods were united to conduct a hydrate exploitation test. 11,12 It is proposed that the combination of depressurization and wellbore heating will greatly increase the production capacity, single thermal excitation is less efficient, the replacement method under high pressure can cause the formation of carbon dioxide hydrate in the wellbore and lead to blockage, and the replacement method under low pressure injection is feasible instead. In 2013, Japan carried out the hydrate exploitation test by depressurization in the condition of 1006 m seawater depth in the sea area of Aichi prefecture, which is the first successful exploration of underwater natural gas hydrate in the world. 13,14 And in 2017, the exploitation test was taken again in the same sea area through the method of depressurization. 15 In 2017, China has carried out a hydrate exploitation test by depressurization in the condition of 1226 m seawater depth in the South China Sea. 16 In general, suitable methods can be used at different stages to exploit the steadier natural gas hydrate reservoir in different reservoir conditions. However, it is not suitable for the metastable and weakly cemented hydrate reservoir in the shallow subsurface of the seafloor, which has no complete trap structure and dense cover. If the above method that separating methane from the hydrate is used to exploit hydrate reservoir in longterm, changes of temperature and pressure field in the reservoir can be caused, then hydrate will decompose in a large amount, reservoir will collapse, and uncontrollable vaporization and dissipation occurs, and then great risks will be brought in turn. [3][4][5] On the basis of this, our research team put forward an innovative technological thinking: the solid fluidization exploitation of marine natural gas hydrate reservoir. The technical process is shown in Figure 1. [3][4][5]17,18 Utilizing the fact that the hydrate can remain stable and decomposes little at low temperature and high pressure in the sea floor, and considering that the metastable and weakly cemented hydrate in shallow subsurface of the seafloor has the characteristics of shallow burial depth and being loose and fragile, the hydrate sediments is exploited in solid excavation method on the seafloor. Then the hydrate sediments are crushed into fine solid particles by secondary crushing, and part of the sand and mud is separated. After that, the particles are mixed with seawater for fluidization, and then are F I G U R E 1 Technical process of solid fluidization exploitation of marine natural gas hydrate reservoir transported in airtight pipeline up to the sea surface platform. Finally, methane gas is obtained by decomposition, separation, and other postprocessing operations.
Obviously, different from the idea in conventional exploitation method that the methane gas is released by hydrates decomposition or replacement in the reservoir, the methane gas is obtained through the controllable and ordered decomposition of hydrate in airtight pipeline and sea surface platform in the process of solid fluidization exploitation of marine natural gas hydrate. This method avoids the original temperature and pressure field change in hydrate reservoir. In this way, a large amount of uncontrollable decomposition and the collapse of the metastable and weakly cemented hydrate reservoir in shallow subsurface of the seafloor are eliminated. It is a kind of innovative method for safe and efficient exploitation of the metastable and weakly cemented hydrate reservoir in shallow subsurface of the seafloor. [3][4][5]17,18 However, the research on the solid fluidization exploitation of marine natural gas hydrate reservoir is in the initial stage, and some theoretical problems need to be further studied. Especially, the temperature in the pipeline (i.e wellbore) is constantly rising and the pressure is constantly decreasing in the process that the hydrate solid particles are transported up to the sea surface platform with the seawater in the airtight pipeline, and the hydrate will decompose at a certain critical position. The large amount of gas produced by hydrate decomposition turns the liquid-solid two-phase flow in wellbore into complex gas-liquid-solid multiphase flow. The hydrate dynamic decomposition particles can be further affected by complex multiphase flow. In this process, the characteristic parameters of the pipe flow are complicated, which is a multiphase non equilibrium pipe flow, and thus the safety risks such as well control and solid particle transportation are increased.
At present, many scholars have carried out the study of multiphase flow regime in wellbore, mainly focusing on flow regime division and prediction, [19][20][21][22] wellbore pressure calculation and measurement. [23][24][25] There is no study on the flow behavior of wellbore multiphase flow with coupling hydrate decomposition. Moreover, the research on hydrate decomposition mainly focuses on decomposition temperature and pressure condition, [26][27][28][29][30] decomposition rate model, and prediction. [31][32][33][34] The dynamic decomposition behavior of hydrate particles in the condition of increasing temperature and decreasing pressure is not studied. Therefore, in this paper we have studied the multiphase non equilibrium pipe flow in the solid fluidization exploitation of marine natural gas hydrate. Mathematical model of wellbore temperature and pressure field, hydrate dynamic decomposition in multiphase riser pipe flow, wellbore multiphase flow coupled hydrate dynamic decomposition are established. The corresponding numerical calculation method is proposed and the accuracy of the models is verified. Then the sensitivity analysis of F I G U R E 2 Physical model for multiphase non equilibrium pipe flow in solid fluidization exploitation of marine natural gas hydrate reservoir | 763 WEI Et al. multiphase non equilibrium flow under different construction parameters is carried out, and the multiphase non equilibrium flow behaviors of solid fluidization exploitation of marine gas hydrate is obtained.

| PHYSICAL MODEL
In solid fluidization exploitation of marine natural gas hydrate, we studied the multiphase non equilibrium pipe flow behaviors in the process of transporting the hydrate solid particles to the sea surface platform with seawater in the airtight pipeline. The simplified physical model is shown in Figure 2. In this process, after the hydrate solid particles entering into bottom hole, they successively past through the subsea horizontal well section, the inclined well section and the vertical well section, and are transported to the sea surface platform finally. In the subsea horizontal well section, the temperature of the mixed fluid in the wellbore is approximately equal to the bottom-seawater temperature, and the wellbore pressure is slightly higher than the bottom-seawater pressure at the same time. Therefore, the hydrate will not decompose in the subsea horizontal well section, and the fluid in the wellbore flows in two phases of liquid and solid. In the inclined well section and the vertical well section, as the well depth decreases, the temperature of the outside seawater increases, the mixed fluid temperature rises in the wellbore, and the wellbore pressure decreases at the same time. Hydrate can only be stable in high pressure and low temperature. Therefore, the hydrate will decompose and a large amount of gas will be released when the hydrate solid particles are transported up to the critical position of decomposition, then the liquid-solid two-phase flow in wellbore changes into complicated gas-liquid-solid multiphase non equilibrium flow.
In this paper, we mainly focus on the study of the behaviors of the wellbore multiphase non equilibrium pipe flow in the rising process of inclined and vertical well. Therefore, we assume that the solid particles in the horizontal segment can be transported safely in order to simplify the model, and the velocity is equal to the velocity of liquid phase. In this study, the mathematical models of wellbore temperature field, pressure field, and hydrate phase equilibrium were established to find out the critical position of the decomposition of hydrate and to determine whether the hydrate solid particles decompose in different well depth. In addition, a mathematical model for hydrate dynamic decomposition in multiphase riser pipe flow was established to calculate the dynamic decomposition rate of hydrate and the rate of gas production. At the same time, a mathematical model for wellbore multiphase flow of coupled hydrate dynamic decomposition was established to calculate the multiphase flow characteristics of the wellbore (gas velocity, liquid velocity, solid velocity, gas holdup, liquid holdup, solid content). In the end, it was necessary to develop a set of numerical calculation method for the multiphase non equilibrium pipe flow in the solid fluidization exploitation, in order to analyze and obtain the relevant behaviors.

| The mathematical model of wellbore temperature field
On the basis of the law of conservation of energy and basic equations of heat transfer, considering the process of solid fluidization exploitation of marine gas hydrate, the temperature distribution model of mixed fluid in wellbore is established: The solution of each parameter is as follows: where T m is the temperature of mixed fluid in the wellbore, K; T w is the temperature of seawater, K; z is the well depth, m; q w is the heat exchange between seawater and wellbore, W m −1 ; q f is the heat generated by the fluid flow friction in the wellbore, W m −1 ; q h is the heat of phase change in the decomposition of hydrate particles during the rising process, W m −1 ; ρ m , ρ g , ρ l, and ρ s are respectively the densities of mixed fluid, gas phase, liquid phase, and solid phase in the wellbore, kg m −3 ; v m , v g , v l, and v s are respectively the velocities of mixed fluid, gas phase, liquid phase and solid phase in the wellbore, m s −1 ; c m , c g , c l, and c s are respectively the specific heat capacities of mixed fluid, gas phase, liquid phase, and solid phase in the wellbore, J kg −1 K −1 ; E g , E l, and E s are respectively the gas holdup, liquid holdup, and solid content, %; D pi and D po are respectively the inner diameter and outer diameter of the wellbore, m; α pi is convection heat transfer coefficient at the inner wall of the wellbore, W m −2 K −1 ; α po is convection heat transfer coefficient at the outer wall of the wellbore, W m −2 K −1 ; λ p is the heat transfer coefficient of wellbore, W m −2 K −1 ; Q m is the volume flow of mixed fluid, m 3 s −1 ; f m is the flow friction coefficient, dimensionless; Δn hyd is the whole amount of substance of the decomposition of hydrate in per length unit and time unit, mol s −1 m −1 ; q hyd is the heat of phase change in per amount of substance hydrate, J mol −1 .
In the Equation (2), seawater temperature Tw can be divided into the mixed layer, the thermocline and the thermostatic layer from top to bottom. 35 According to the actual situation in South China Sea, 36

| The mathematical model of wellbore pressure field
The wellbore pressure distribution model is established based on the process of multiphase mixed fluid flow in wellbore in solid fluidization exploitation of marine gas hydrate: The solution of each parameter is as follows: where p m is the wellbore pressure, MPa; p G , p F, and p A are respectively the pressure drops resulting from gravity, frictional resistance and a change in fluid velocity, MPa m −1 ; g is the acceleration of gravity, m s −2 ; θ is the hole drift angle, radian.

| The phase equilibrium model of hydrate
On the basis of the process of solid fluidization exploitation of marine gas hydrate, a hydrate phase equilibrium model is required to be established in order to determine whether the hydrate solid particle decomposes in a certain temperature and pressure condition during the rising process. In order to simplify the calculation, this paper adopts the hydrate phase equilibrium model established through experimental data analysis by Dzyuba et al. 37 where p eq is the hydrate phase equilibrium pressure at a certain temperature of T m , MPa.
According to Equation (10), the equilibrium pressure p eq of hydrate can be calculated in view of the wellbore temperature of T m when the hydrate solid particles rise to a certain position in the wellbore. If p eq is less than the wellbore pressure p m at this time, the hydrate will not decompose. If p eq is equal to the wellbore pressure p m at this time, it is regarded as the critical position of hydrate decomposition. If p eq is greater than the wellbore pressure p m at this time, the hydrate will decompose.

| The mathematical model of hydrate dynamic decomposition in multiphase riser pipe flow
In the riser pipe flow, according to the equilibrium model of hydrate, the hydrate solid particles will not decompose before rising to the critical decomposition position, and the decomposition rate is zero in this process. The hydrate solid particles will decompose after arriving to the critical decomposition position. In order to calculate the decomposition of hydrate, based on Kim model, 31 the dynamic decomposition model of hydrate in multiphase riser pipe flow was established, in which the gas hydrate is assumed to be methane hydrate.
where n hyd is the amount of hydrate substances in solid particles, mol; t hyd is the hydrate decomposition time, s; k hyd is the hydrate decomposition rate constant, mol s −1 m −2 MPa −1 ; S hyd is the decomposition surface area of hydrate solid .6339 ln p eq + 264.9661 particle, m 2 ; f eq and f m are respectively the fugacities of methane gas at temperature of T m and pressure of p eq and p m , MPa. The solution of each parameter is as follows.

| The hydrate decomposition rate constant
The hydrate decomposition rate constant has the following relation: where k hydc is the decomposition reaction rate of hydrate itself, mol s −1 m −2 MPa −1 ; k hydf is the mass transfer rate of methane gas, mol s −1 m −2 MPa −1 .
According to the Arrhenius equation, 38 the decomposition reaction rate of hydrate itself is calculated as follows: where k hydc 0 is the intrinsic decomposition rate constant of hydrate, and it is valued as 8.77 × 10 11 mol s −1 m −2 MPa −1 according to Kim 31 ; E act is the activation energy of hydrate decomposition, and it is valued as 7.83 × 10 4 J mol −1 according to Kim 31 ; R is the general gas constant, valued as 8.314 J mol −1 K −1 .
In solid fluidization exploitation of marine natural gas hydrate, the hydrate sediments were crushed into fine solid particles, and it is assumed that the hydrate solid particles are homogeneous spheres. In the riser pipe flow, the hydrate solid particles are transported up. The hydrate decomposition rate constant is influenced by the decomposition reaction rate of hydrate itself and the mass transfer rate of methane gas. The mass transfer of methane gas is from the hydrate solid particles to the pipe flow. So the mass transfer from the sphere to the fluid can be used in this process, and they are valid in this model. The empirical equation to describe the mass transfer is written as 39 : where n sh , Re f, and n sc are respectively Sherwood number, Reynolds number, and Schmidt number of mass transfer. According to the theory of mass transfer theory, 40 they can be written as Thus, the mass transfer rate of methane gas is obtained: where d s is the diameter of hydrate solid particle, m; D AB is the diffusion coefficient of methane gas to pipe flow, m 2 s −1 ; μ m is the viscosity of the mixed fluid in the wellbore, Pa s; Z g is the compression factor of natural gas, dimensionless.

| The decomposition surface area of hydrate solid particle
It is assumed that the decomposition of hydrate solid particles occurs at the interface only. 41 The decomposition surface area of hydrate solid particle is written as and where d hyd is the equivalent diameter of hydrate solid particles, m; φ is sphericity, dimensionless; V hyd is the volume of hydrate in solid particle, m 3 ; M hyd is the molar mass of hydrate, kg mol −1 ; ρ hyd is the density of the hydrate, kg m −3 ; E hyd is the abundance of hydrate in solid particle, %; V s is the volume of hydrate solid particle, m 3 .

| The fugacity of methane gas
According to the definition of gas fugacity, 42 it can be concluded that: where R atm is the general gas constant expressed in standard atmosphere, valued as 0.08206 atm L mol −1 K −1 ; T x is the environmental temperature, K; p x is environmental pressure, atm; p* is the reference pressure, atm; f x is the fugacity of methane gas at T x and p x , atm; V x is the molar volume of methane gas at T x and p x , L mol −1 .
Suppose that the methane gas state is subject to the R-K state equation, 39 the molar volume V x of methane gas under temperature of T x and pressure of p x has the following relation: where a and b are R-K constants, a provides a measure of intermolecular attraction, b provides the size of the molecule. The calculation is as follows 39 : where T c is the critical temperature of methane gas, K; p c is the critical pressure of methane gas, atm.
Therefore, the calculation equation of methane gas fugacity is derived: In addition, the fugacity of methane gas can be obtained at temperature of T m and pressure of p eq and p m.
where V m and V eq are the molar volumes of methane gas at temperature of T m and pressure of p eq and p m respectively, L mol -1 .

| The mathematical model of wellbore multiphase flow coupled hydrate dynamic decomposition
The liquid and solid phases were considered as incompressible flow in the axial and radial directions. Therefore, the wellbore multiphase flow models coupled hydrate dynamic decomposition is established as follows.
Mass conservation equation of gas and liquid phases is written as: Mixed momentum equation of gas, liquid, and solid phases in wellbore is written as: where m hydg and m hydl are respectively the mass change rates of gas and liquid phases in unit length, kg s −1 m −1 . In the riser pipe flow, the hydrate solid particles will not decompose before rising to the decomposition critical position, the decomposition rate is zero in this process, and m hydg and m hydl are both zero. The hydrate solid particles will decompose after arriving to the critical position. Then m hydg and m hydl are calculated as follows: where Δn hyd is the total amount of substance of hydrate decomposition in per unit time and unit length, mol s −1 m −1 ; V methane and V water are respectively the volumes of gas and liquid production of 1 unit volume hydrate after decomposition in the standard conditions, and reference has been made to take 164 and 1 unit volume, respectively 43 ; ρ methane and ρ water are the densities of methane gas and water in standard conditions, respectively, kg m −3 .
In the wellbore multiphase flow models, the liquid velocity depends largely on the liquid delivery capacity, and the velocity equations of gas and solid are written as 44,45 : where C g is distribution parameter of gas velocity, dimensionless; v gt is slip velocity of gas phase, which is related to flow regime, m s −1 ; C s is distribution parameter of solid velocity, dimensionless; v st is settling velocity of solid phase, which is related to flow regime, m s −1 . C g , v gt , C s, and v st in this study are calculated by the method in Reference. 45 In addition, where and v sg is the apparent gas velocity, m s −1 ; v sl is the apparent liquid velocity, m s −1 ; σ gl is the interfacial tension between gas and liquid, which is calculated by the equation based on experimental data, 50 N m −1 ; v 0∞ is the rising velocity of a single bubble in an infinite medium, m s −1 . It can be seen from Equations 30-36 that the flow regime partition are mainly affected by apparent gas velocity, apparent liquid velocity and mixed fluid velocity. In addition, the value of solid content directly affects the mixed fluid velocity in wellbore. Therefore, the influence of solid content on flow regime partition is considered through the change of mixed fluid velocity.

| The auxiliary mathematical model
In the mathematical model of wellbore temperature field, the general expression of convection heat transfer coefficient at the inner and outer walls of the wellbore are defined as: where Nu pi and Nu po are Nusselt numbers at the inner and outer walls of the wellbore, dimensionless. The solution of them are as follows 51-53 : where Re pi and Re po are Reynolds numbers at the inner and outer walls of the wellbore, dimensionless; λ m and λ w are the heat transfer coefficients of the mixed fluid in the wellbore and the seawater outside the wellbore, W m −2 K −1 ; c w , v w , ρ w, and μ w are respectively the specific heat capacity, velocity, density, and viscosity of the seawater outside the wellbore, and the physical units are respectively J kg −1 K −1 , m s −1 , kg m −3 , and Pa s; a x and b x are calculating coefficients that is identified by Reynolds number, dimensionless. 52 The viscosity of the mixed fluid in the wellbore is written as where μ g and μ l are gas and liquid viscosities, respectively, Pa s. The gas viscosity is calculated by the equation based on a large amount of experimental data by Lee et al. 54 : where δ g is the gas relative density to water, dimensionless. The liquid viscosity is calculated by the equation based on experimental data by Beggs and Brill 50 : The flow friction coefficient is dependent on the flow state, which is identified by Reynolds number 55 : and where ξ is the roughness of the inner wall of the wellbore, m.

| Boundary and initial conditions
In order to study the multiphase non equilibrium flow behaviors in the solid fluidization exploitation of marine natural gas hydrate, the finite difference numerical calculation method is used in this paper. Meanwhile, the influential factors of the boundary and initial conditions among each control region are taken into account in the solving model.
In the process of solid fluidization exploitation of marine natural gas hydrate reservoir, the temperature of mixed fluid at bottom hole is equal to the injection seawater temperature at the seafloor. Therefore, the temperature of the mixed fluid in wellbore in the subsea horizontal well section remains the same as the seawater temperature at the outside seafloor: The pressure at the wellhead is equal to the wellhead back pressure: where z H is the vertical depth of the subsea horizontal well section, m; T m (z H ,t) and T w (z H ) are respectively the temperatures of mixed fluid in the wellbore of subsea horizontal section and the external seawater, K; p m (0,t) and p b are respectively the pressures at the wellhead and the wellhead back pressure, MPa.
In the numerical calculations of this study, in the well section that before hydrate decomposition, we consider the hydrate as a single phase, the solid particle in the wellbore is composed of the adglutinate of hydrate and silt, and the density of solid particle is equal to the volumetrically weighted average of their densities. In the well section during hydrate decomposition, the hydrate decomposes into gas and water, the solid particle in the wellbore is composed of residual silt and the adglutinate of hydrate and silt, and the density of solid particle is equal to the volumetrically weighted average of their densities. In the well section after hydrate decomposition, the hydrate has completely decomposed into gas and water, the solid particle in the wellbore is composed of residual silt, and the density of solid particle is equal to that of residual silt.

| Numerical calculation method
In the process of finite difference numerical calculation, the spatial domain is the wellbore, and the time domain is the hydrate solid particles transported from the bottom hole to the wellhead, and the calculation sequence is from the wellhead to the bottom hole.
Suppose at time n, the parameters of any two nodes i and i + 1 in the wellbore are known conditions. The numerical calculation process is illustrated with the multiphase non equilibrium pipe flow process from n to n + 1 of node i and i + 1 (shown in Figure 3). α, β, γ in Figure 3 are all calculated error accuracy.

| Verification of numerical calculation
In order to verify the accuracy of the numerical calculation method established in this paper, the large physical simulation experimental results of marine natural gas hydrate solid fluidization exploitation by Zhou et al. 5,43 was used to be compared. The function of the experimental facility of solid fluidization exploitation is to simulate the whole technological process. The experimental facility has a vertical plexiglass pipe of 30 m and a horizontal plexiglass pipe of 56 m. In addition, through the continuous pressure and temperature regulation, it can divide the process of upward flow into several 30 m-length simulations. Finally, the simulation of the whole process can be completed by combining the results of each 30 m-length simulation. More introductions of this experimental facility are in the Reference. 43 The basic experimental data in Zhou's paper are as follows: the depth of the seawater is 1315 m, the well depth is 1500 m, the diameter of the experimental tube is 0.076 m, the simulated surface temperature of sea is 308 K, the volume fraction of hydrate in solid particle is 70%, the simulated delivery capacity is 0.008 m 3 s −1 , the density of liquid phase is 1030 kg m −3 , the solid throughput is 0.00015 m 3 s −1 . The values of the roughness of the inner wall of the wellbore is set to be 0.00001 m because of the plexiglass pipe, and the values of other related parameters refer to Table 1 and some are calculated by the mathematical models.
The numerical simulation method in this paper is applied to compare with the experimental data in Zhou's paper. The wellbore temperature, pressure, gas holdup, liquid holdup, solid content, gas velocity, liquid velocity, and solid velocity curve are obtained as shown in Figures 4 and 5.
In Figures 4 and 5, compared with the experimental results in Zhou's paper, we can see that: the error between the numerical results obtained in this paper and the experimental results in Zhou's paper is small, and the variation trend is consistent. From the analysis in the physical model of solid fluidization exploitation of marine natural gas hydrate, we can see that the wellbore flow from bottom to wellhead can be divided into the liquid-solid two-phase flow before hydrate decomposition, the gas-liquid-solid multiphase flow during hydrate decomposition, and the gas-liquid-solid multiphase flow after hydrate decomposition. It is obvious that the characteristic parameters of the pipe flow during hydrate decomposition are the most complex, which is a multiphase non equilibrium pipe flow. Furthermore, the numerical calculations for this process are the most complex. As can be seen from Figure 5a, the gas holdup gradually increases from 0, which indicates that the hydrate just begins to decompose in the simulating well section from 670 m to 640 m. At this point, the wellbore flow is gas-liquid-solid multiphase flow during hydrate decomposition. Then, by comparing the numerical and experimental results of the most complex well section in the whole pipe flow process, the accuracy of the whole set of numerical calculation method can be verified better, and that will be more reasonable in the application of other well sections. Therefore, the accuracy of the numerical model and the calculation method established in this paper has been verified, which lays the foundation for the application of the numerical model below.

APPLICATIONS
In the process of solid fluidization exploitation of marine natural gas hydrate reservoir, the well structure is shown in Figure 2: the well depth is 1400 m, the vertical depth is 1000 m, the vertical well section is 0~830 m, the inclined well section is 830~1100 m, the horizontal well section is 1100~1400 m, the outer diameter of wellbore is 0.508 m, and the inner diameter is 0.476 m. The characteristics of natural gas hydrate reservoir in the South China Sea are 3-5 : the seawater depth is 1000 m, the sea surface temperature is 293 K, the seawater density is 1030 kg m −3 , the volume fraction of hydrate in solid particle is 70%, the hydrate density is 910 kg m −3 , the silt density in the hydrate sediment is 2500 kg m −3 , the hydrate solid particle density is the weighted average of hydrate and silt densities and it is 1387 kg m −3 , the average sphericity of hydrate solid particle is 0.8, and the average diameter is 0.005 m.
In the numerical calculation, the value of gas density is closely related to temperature and pressure, which can be obtained by the equation of gas state. We assume that the liquid and solid densities do not change with temperature and pressure. The calculations of gas holdup, liquid holdup and solid content depend on the flow regime, and the calculation processes refer to the Reference. 50 In addition, the values of other related parameters are listed in Table 1. [44][45][46][47][48][49][50][51][52][53][54][55] In order to obtain the multiphase non equilibrium pipe flow behaviors in the process of transporting the hydrate solid particles to the sea surface platform with seawater in the airtight pipeline, the established mathematical model, and numerical calculation method are used to carry out the numerical model application analysis of different liquid delivery capacities, solid throughputs, liquid densities, and wellhead back pressures.

| Analysis under different liquid delivery capacities
The numerical calculation parameters are as follows: the liquid delivery capacities are 0.03 m 3 s −1 , 0.05 m 3 s −1 , and 0.07 m 3 s −1 , the solid throughput is 1.5 m 3 h −1 , the liquid density is 1030 kg m −3 , and the wellhead back pressure is 0.1 MPa. Through numerical calculation under different liquid delivery capacities in the solid fluidization exploitation of marine natural gas hydrate reservoir, the wellbore temperature and pressure are obtained in Figure 6, hydrate 770 | WEI Et al. dynamic decomposition behaviors in multiphase riser pipe flow in Figure 7, and wellbore multiphase flow behaviors of coupled hydrate dynamic decomposition are obtained in Figure 8.
As can be seen from Figure 6, in the process of transporting the hydrate solid particles to the sea surface platform with seawater, the wellbore temperature at the subsea horizontal well section from 1400 m to 1100 m is equal to the temperature of outside seawater. Then at the inclined and vertical well section from 1100 m to 0 m, the wellbore temperature T m is increased by the increase of seawater temperature T w .
In addition, the heat exchange time of mixed fluid in the wellbore and seawater are shortened with the increase of liquid delivery capacity (from 0.03 to 0.07 m 3 s −1 ), and the heat transferred from seawater to the wellbore decreases. Therefore, the temperature of the wellbore is reduced, and then the phase equilibrium pressure p eq calculated by the hydrate phase equilibrium model (Equation (10)) is reduced (Figure 6b).
With the increase in liquid delivery capacity, the solid content E s in the wellbore decreases under the same solid throughput (Figure 8f), resulting in a decrease in the mixed  fluid density ρ m , which tends to decrease the wellbore pressure p m . However, the increase in each phase velocities (Figure 8a-c) will increase the friction pressure drop p F and tend to increase the wellbore pressure p m . Therefore, the change in wellbore pressure p m is not obvious in Figure 6b. In Figure 6b, the intersection point of the wellbore pressure curve and the phase equilibrium pressure curve in the condition of any definite liquid delivery capacity is the critical position of the hydrate solid particles decomposition in the upward pipe flow, corresponding to Figure 7 is the position where the amount of substance n hyd of the hydrate in solid particles starts to decrease. After that, the amount of substance decreases as the hydrate solid particles rise continuously in the wellbore in Figure 7. The position where the amount of substance comes at 0 is where the hydrate completely decomposes. At the same time, it can be seen that, with the increase of the delivery capacity of liquid phase (from 0.03 to 0.07 m 3 s −1 ), the critical position of hydrate decomposition moves up, and the position of complete decomposition also moves up.
According to the assumption of the physical model, at the subsea horizontal well section from 1400 m to 1100 m, the velocity of solid particles is equal to the velocity of liquid phase, and the fluid in the wellbore flows in liquid-solid two phases because of no decomposition of hydrate. Thus, as can be seen from Figure 8a-c, the liquid velocity v l and the solid velocity v s of any definite liquid delivery capacity are equal, and the gas phase velocity v g is 0 at the subsea horizontal well section. Then at the inclined well section from 1100 m to 830 m, the liquid velocity v l remains unchanged because the section area of the wellbore is unchanged (the inner diameter of wellbore is constant), and the solid velocity v s is gradually decreased under the negative influence of gravity. At the vertical well section from 830 m to 0, in combination with the critical position of hydrate decomposition in Figures 6b and 7, the gas velocity is still 0 when the hydrate does not decompose, and the liquid velocity v l and solid velocity v s remain stable. When the hydrate solid particles are transported up to the critical position, the hydrate in the solid particles decomposes into gas and water. Therefore, the twophase flow of liquid and solid in wellbore becomes complex gas-liquid-solid multiphase flow. As the multiphase flow continues to rise, the well depth decreases, the gas expands and the volume increases, and the gas velocity v g increases (Figure 8a), and then its carrying effect increases the liquid velocity v l (Figure 8b). However, due to the decrease in the density of mixed fluid in the wellbore, the carrying capacity decreases, so the solid velocity v s decreases with the well depth (Figure 8c).
At the same time, when compared with the effect of different liquid delivery capacities in Figure 8a-c, we can see that the velocities of gas, liquid, and solid in wellbore increase with the increase of liquid delivery capacity (from 0.03 m 3 s −1 to 0.07 m 3 s −1 ).
In Figure 8d-f, combining 8a-c we can see that at the subsea horizontal well section from 1400 m to 1100 m, the gas holdup E g is 0, the liquid holdup E l , and solid content E s remain stable. Then at the inclined well section from 1100 m to 830 m, the solid velocity is gradually decreased, causing the solid content E s increase gradually, and the corresponding liquid holdup E l in the wellbore is decreased. At the vertical well section from 830 m to 0 m, combining Figure 7, when the hydrate does not decompose, the gas holdup E g is 0, and the liquid holdup E l and the solid content E s remain stable. When the hydrate solid particles are transported up to the critical decomposition position, the hydrate in the solid particles decomposes into gas and water. Therefore, the gas holdup E g increases, the liquid holdup E l increases, and the solid content E s decreases in the process of the hydrate decomposition. The liquid holdup E l does not increase anymore and the solid content E s does not decrease anymore when the hydrate solid particles are transported up to the full decomposition position. In addition, the gas expansion volume increases and the gas holdup increases as the well depth decreases, resulting in a decrease in the density of mixed fluid in the wellbore, and a decrease in the carrying capacity. Therefore, the solid content E s increases, and the corresponding holdup E l decreases.
At the same time, while comparing the effect of different liquid delivery capacities in Figure 8d-f, we can see that with the increase in liquid delivery capacity (from 0.03 m 3 s −1 to 0.07 m 3 s −1 ), the gas holdup in the wellbore decreases, the liquid holdup increases, and the solid content decreases. When the liquid delivery capacity comes to its minimum as 0.03 m 3 s −1 , the solid velocity v s suddenly rises (Figure 8c), and the solid content E s suddenly decreases (Figure 8f) near the wellhead (the well depth is close to 0). This shows that the high gas holdup at the wellhead causes the transformation of the flow regime of multiphase flow when other conditions Combining Figures 6-8, we can see that when the liquid delivery capacity is lower, the transportation capacity for solid particles is lower, the solid content is higher, and the risk of plugging system is higher, the gas velocity and gas holdup are higher, the flow regime transformation of multiphase flow occurs more easily, the solid velocity increases more and the solid content decreases more near the wellhead, and the well control safety risk is more serious. So the liquid delivery capacity should be increased appropriately in order to prevent well control problems and improve the transportation capacity for solid particles in the field construction of solid fluidization exploitation of marine natural gas hydrate reservoir.

| Analysis under different solid throughputs
In the hydrate exploitation test in the South China Sea by China, the production of natural gas lasted 60 days, the cumulative production is 309 000 m 3 , and the average daily production is 5150 m 3 d −1 . 16 According to this engineering background, we set the solid throughputs in numerical calculation to be 1.  Figure 9, hydrate dynamic decomposition behaviors in multiphase riser pipe flow in Figure 10, and wellbore multiphase flow behaviors of coupled hydrate dynamic decomposition are obtained in Figure 11.
As can be seen from Figures 9-11, similar to the process of . Therefore, we only analyze the effect laws of different solid throughputs in this section.
In Figure 9, the amount of the heat transferred from seawater to the wellbore does not change much with the increase in solid throughput (from 1.5 m 3 h −1 to 4.5 m 3 h −1 ). Therefore, it is not obvious that the wellbore temperature T m changes, and the calculated phase equilibrium pressure p eq changes little.
In Figure 9b, the higher solid throughput will bring more solid phase into the wellbore in per time unit, and will create a rise trend of the density of mixed fluid in the wellbore. The fluid in the wellbore flows in liquid-solid two phases in the lower well section before the hydrate decomposes. Therefore, the density of the mixed fluid in the wellbore increases, resulting in the increase of the gravity pressure drop p G , thus increasing the wellbore pressure p m in the lower well section. However, the gas phase is released from the hydrate in the upper well section where the hydrate has decomposed, and the wellbore flow becomes in gas-liquid-solid multiphase. The gas phase has a tendency to decrease the density of mixed fluid in the wellbore, resulting in the decrease in the gravity pressure drop p G . Therefore, the wellbore pressure p m does not change obviously with the increase in solid throughput in the upper well section in combination of the effect of the solid phase. At the same time, the critical decomposition position of the hydrate solid particles in the upward transportation with the seawater (the intersection of the curve of wellbore pressure p m and the curve of phase equilibrium pressure p eq ) does not change obviously with the increase in solid throughput, and the complete decomposition position does not change obviously too ( Figure 10).
In Figure 11, with the increase in solid throughput (from 1.5 m 3 h −1 to 4.5 m 3 h −1 ), more solids enter into the wellbore in per time unit, and the higher solid content is (Figure 11f). When the solid particles are transported upward to the critical decomposition position with seawater, the higher the amount of gas produced by the hydrate decomposition is, the higher the gas holdup E g (Figure 11d) will be, and the lower the corresponding liquid holdup E l will be. With the decrease in well depth, the larger the solid throughput is, the larger the volume of gas after expansion will be, the higher the gas velocity v g will be, and the liquid velocity v l is higher (Figure 11b) by its carrying effect. Meanwhile, the gas holdup at the wellhead is higher in the condition of large solid throughput, resulting in a more severe transformation from the bubble flow or dispersed bubble flow in the lower well section to the slug flow, which leads to the greater transformation of the carrying capacity. Therefore, the more solid velocity increases at the wellhead (Figure 11c), the more solid content decreases (Figure 11f). And it poses a greater challenge to well control safety at the same time.
As can be seen from the combination of Figures 9-11, increasing the solid throughput can increase the amount of hydrate transported to the sea surface platform in per time unit, and further increase the natural gas production in the field construction of solid fluidization exploitation of marine natural gas hydrate reservoir. However, the solid content will increase and the risk of plugging system will increase. Then the amount of gas produced by the hydrate decomposition will be higher and the gas holdup will be higher, the flow regime transformation of multiphase flow will occur more easily, the solid velocity will increase more and the solid content will decrease more near the wellhead, and the well control safety risk will be more serious. Therefore, for avoiding the risks of plugging system and well control safety, the solid throughput can only be improved appropriately.

| Analysis under different liquid densities
The numerical calculation parameters are as follows: the liquid densities are 1030 kg m −3 , 1130 kg m −3 and 1230 kg m −3 , the liquid delivery capacity is 0.03 m 3 s −1 , the solid throughput is 1.5 m 3 h −1 , and the wellhead back pressure is 0.1 MPa. Through numerical calculation under different liquid densities in the solid fluidization exploitation of marine natural gas hydrate reservoir, the wellbore temperature and pressure are obtained in Figure 12, hydrate dynamic decomposition behaviors in multiphase riser pipe flow in Figure 13, and wellbore multiphase flow behaviors of coupled hydrate dynamic decomposition are obtained in Figure 14.
The change laws of each parameter in Figures 12-14 from the bottom hole to the wellhead are similar to the analysis in the section of "Application analysis under different liquid delivery capacities." Therefore, we analyze the effect laws only of different liquid densities in this section.
In Figures 12 and 13, the heat transferred from seawater to wellbore is not significant with the increase in liquid density (from 1030 kg m −3 to 1230 kg m −3 ). Therefore, the change in wellbore temperature T m is not obvious, and the calculated phase equilibrium pressure p eq is not changed much. In the wellbore, the density of mixed fluid increases, which increases the gravity pressure drop p G greatly, and then increases the wellbore pressure p m . The critical decomposition position (the intersection of the curve of wellbore pressure p m and the curve of phase equilibrium pressure p eq ) and the complete decomposition position of the hydrate solid particles in the upward transportation with the seawater both move up.
In Figure 14, the velocity of hydrate solid particles in the subsea horizontal section from 1400 m to 1100 m is equal to the velocity of liquid phase according to the assumption of the physical model. Therefore, both the liquid velocity v l and the solid velocity v s are not changed with the increase in liquid density (Figure 14b,c). In the inclined and vertical well sections from 1100 m to 0 m, the higher the liquid density is, the higher the solid carrying capability of the wellbore fluid, the greater the solid velocity v s (Figure 14c), the lower the solid content E s (Figure 14f), and the higher the liquid holdup E l (Figure 14e) will be, and liquid velocity v l will stay unchanged ( Figure 14b). Meanwhile, the wellbore pressure p m is higher (Figure 12b) in the condition of high liquid density, which will cause less expansion of the gas after the hydrate solid particles are transported upward to the critical decomposition position. Therefore, the gas velocity v g and the gas holdup E g are slightly lower (Figure 14a,d). At the same time, it is showed that the transformation of the flow regime from the bubble flow or dispersed bubble flow to the slug flow is more difficult in the higher liquid density condition, the solid velocity v s and the solid content E s (Figure 14c,f) at the wellhead are changed less, and it is more helpful to ensure the well control safety.
Combining Figures 12-14, we can see that increasing the liquid density can significantly improve the transportation capacity for solid hydrate particles, and reduce the solid content and the risk of plugging system. It can also reduce the gas velocity and gas holdup, and make the flow regime of multiphase flow steadier. Then the solid velocity and solid content change smaller near the wellhead, and the well control safety risk is lower. Therefore, it is necessary to increase the liquid density appropriately in the field construction of the solid fluidization exploitation of marine natural gas hydrate reservoir.

| Analysis under different wellhead back pressures
The numerical calculation parameters are as follows: the wellhead back pressures are 0.1 MPa, 1 MPa, and 2 MPa, the liquid delivery capacity is 0.03 m 3 s −1 , the solid throughput is 1.5 m 3 h −1 , and the liquid density is 1030 kg m −3 . Through numerical calculation under different wellhead back pressures in the solid fluidization exploitation of marine natural gas hydrate reservoir, the wellbore temperature and pressure are obtained in Figure 15, hydrate dynamic decomposition behaviors in multiphase riser pipe flow in Figure 16, and wellbore multiphase flow behaviors of coupled hydrate dynamic decomposition are obtained in Figure 17.
The change laws of each parameter in Figures 15-17 from the bottom hole to the wellhead are similar to the analysis in the section of "Application analysis under different liquid delivery capacities". Therefore, we analyze the effect laws only of different wellhead back pressures in this section.
In Figures 15 and 16, the heat transferred from seawater to wellbore is not significant with the increase in wellhead back pressure (from 0.1 MPa to 2 MPa). Therefore, the change in wellbore temperature T m is not obvious, and the calculated phase equilibrium pressure p eq is not changed much. The In Figure 17, before the hydrate solid particles are transported upward to the critical decomposition position with seawater, the fluid in the wellbore flows in liquid-solid two phases, the velocity of each phase (v g , v l , and v s ), and the content of each phase (E g , E l , and E s ) all does not unchanged with the increase in wellhead pressure (from 0.1 MPa to 2 MPa). Then after the hydrate solid particles are transported upward to the critical decomposition position, the hydrate decomposes to produce gas. The higher the wellhead back pressure is, the smaller the gas expansion, the lower the gas velocity v g and the gas holdup E g (Figure 17a,d), the lower the liquid velocity v l (Figure 17b), the higher the liquid holdup E l (Figure 17e) will be. Meanwhile, the higher the density of mixed fluid in the wellbore is, the greater the solid carrying capability, the larger the solid velocity v s (Figure 17c), and the lower the solid content E s (Figure 17f) will be. At the same time, it is showed that the transformation of the flow regime from the bubble flow or dispersed bubble flow to the slug flow is more difficult in the higher wellhead back pressure condition, the solid velocity v s and the solid content E s (Figure 17c,f) at the wellhead are changed less, and it is more helpful to ensure the well control safety.
Combining Figures 15-17, we can see that increasing the wellhead back pressure can significantly improve the wellbore pressure. Then the gas expansion will be smaller, the gas velocity and gas holdup will be lower, the flow regime of multiphase flow will be steadier, and the well control safety risk will be lower. Therefore, it is necessary to increase the wellhead back pressure appropriately in the field construction of the solid fluidization exploitation of marine natural gas hydrate reservoir.

| CONCLUSIONS
In this paper, the physical model and the system mathematical models of the multiphase non equilibrium pipe flow in the solid fluidization exploitation of marine natural gas hydrate reservoir are established: mathematical model of wellbore temperature field and pressure field, mathematical model of the phase equilibrium of hydrate, mathematical model of the dynamic decomposition of natural gas hydrate in the multiphase riser pipe flow, and the mathematical model of wellbore multiphase flow coupled hydrate dynamic decomposition. Meanwhile, the numerical calculation methods for the multiphase non equilibrium pipe flow in the solid fluidization exploitation of marine natural gas hydrate reservoir is developed and verified. In addition, the numerical model is applied to analyze and obtain the flow behaviors of the multiphase non equilibrium flow in the solid fluidization exploitation of marine natural gas hydrate reservoir: 1. The increase in liquid delivery capacity can significantly reduce the wellbore temperature and phase equilibrium pressure of the hydrate, so that the critical decomposition position of hydrate moves up. It can also significantly improve the velocities of gas, liquid, and solid phase, and can reduce the gas holdup and solid phase content. 2. With the increase in solid throughput, the wellbore temperature and hydrate phase equilibrium pressure do not change, the wellbore pressure increases in the lower well section and does not change obviously in the upper well section, the critical decomposition position of hydrate stays constant, the solid content and gas holdup at wellhead increase significantly, the liquid holdup decreases significantly, and the velocities of gas, liquid, and solid phases all increase. 3. The increase in liquid density can significantly increase the wellbore pressure, so that the critical decomposition position of hydrate moves up. It can also significantly improve the solid velocity, so that the solid content is reduced, and the gas velocity and gas holdup decrease slightly and the liquid holdup increases at the same time. 4. The increase in wellhead back pressure can significantly increase the wellbore pressure, so that the hydrate critical decomposition position moves up. It can also significantly reduce the gas holdup after the hydrate decomposes, reduce the solid velocity and liquid velocity, increase the solid velocity, reduce the solid content, and increase the liquid holdup. 5. The appropriate increase in solid throughput can increase the production of natural gas in the field construction of the solid fluidization exploitation of marine natural gas hydrate reservoir. However, problems such as well control risks will intensify. In order to prevent wellbore flow safety problems such as well control risk, and improve the transportation capacity for hydrate solid particles, the liquid delivery capacity, liquid density, and wellhead back pressure should be increased appropriately.