DOI : https://doi.org/10.5281/zenodo.18889692
- Open Access
- Authors : Chinedu Charles Alamezie
- Paper ID : IJERTV15IS020795
- Volume & Issue : Volume 15, Issue 02 , February – 2026
- Published (First Online): 07-03-2026
- ISSN (Online) : 2278-0181
- Publisher Name : IJERT
- License:
This work is licensed under a Creative Commons Attribution 4.0 International License
Comparing Analytical and Numerical Solutions to A Comprehensive Hydrate Kinetic Growth Model
Chinedu Charles Alamezie University of Aberdeen
Abstract – Gas hydrate formation in subsea pipelines poses a significant risk to flow assurance, with deposits capable of restricting or completely blocking hydrocarbon transport. While the thermodynamic stability regions of methane hydrates are well established, accurately predicting their kinetic growth remains challenging due to a limited knowledge and research around the controlling mechanisms. Intrinsic kinetics dominate the early stages of hydrate formation, with heat and mass transfer constraints becoming increasingly important as growth progresses. Existing kinetic growth models do not account for the summation of all influencing growth mechanism hence necessitating this work. To address this issue, we employ the transport equation to complete an already existing growth model. The second-order partial differential nature of the new models expression together with its practical inhomogeneous boundary conditions needs a solution. This new model is then solved via two solution techniques: analytical and numerical. Firstly, in the analytical solution, we employ the separable variable technique to solve our model, while for the numerical solution, Taylors series expansion and the Crank-Nicolson method were implemented for increased stability and accuracy. The numerical and analytical solutions to the proposed comprehensive model are then implemented and compared, with discrepancies explained. The work delivers novel growth expressions and enhances understanding of the coupled effects of intrinsic kinetics, heat transfer, and mass transfer in methane hydrate formation. The outcomes provide both theoretical insight and practical tools to support the design of hydrate management strategies in production systems.
Keywords: gas hydrate growth, Crank-Nicolson method, separation of variables, inhomogenous boundary condition, growth mechanisms
INTRODUCTION
The shift into deeper waters requires longer subsea tiebacks to transport the hydrocarbon fluids from the wellhead to surface processing facilities. In these transportation flowlines, heat exchange with the environment, Joule-Thompson cooling, or a combination of both cause the subsea high-pressure flowlines to cool the fluids to within the thermodynamic stability region of hydrates. Natural gas hydrate formation occurs when associated water and small natural gas molecules, typically methane, ethane, ethylene, iso-propane, and carbon dioxide, are present at high-pressure, low-temperature conditions within the flowlines. In the non-stoichiometric hydrate crystal lattice, the water molecules form hydrogen-bonded cage-like structures that trap the gas (Creek et al., 2011; Kelkar et al., 1998; Max et al., 2006; Stoporev et al., 2018). The formation of gas hydrates causes increased pressure- drop along pipe sections and will ultimately plug the flowline if not addressed. Blockage of flowlines may raise safety and integrity concerns, as well as result in financial losses due to lost production time and remedial operations. It is, therefore, imperative to control hydrate formation. There is a growing interest in managing the build-up of hydrate plugs rather than preventing their appearance due to the high cost of preventive measures.
Gas hydrate formation is considered a two-stage process: nucleation and growth. The hydrate nucleation process involves clustering of water molecules around molecules of gas to create crystal embryos whose stability depends on local mass, temperature, and pressure fluctuations (Kashchiev & Firoozabadi, 2002a; Ke et al., 2019). Hydrate growth is the process in which self-sustainable hydrate nuclei continue to grow. The boundary between nucleation and growth is challenging to place, as they both occur at a molecular level (Ke et al., 2019; Stoporev et al., 2018). Experimental studies detect the onset of hydrates macroscopically, further observing the ensuing hydrate growth. The study of natural gas hydrate nucleation processes in deepwater gas pipeline system conditions is new and increasingly gaining relevance industrially and scientifically (Arjun & Bolhuis, 2020; Ke et al., 2019).
Three controlling mechanisms is generally accepted to govern hydrate growth: The intrinsic growth kinetics mechanism, which relies on the reaction kinetics at the hydrate surface, the heat transfer mechanism, which is dependent on the exothermic nature of hydrate formation and convective heat dissipation, and the mass transfer mechanism, which involves the mass transfer of water and guest molecules to the growing surface (Yin et al., 2016; Yin et al., 2018; Ke et al., 2019; Filarsky et al., 2023). Pioneering work on hydrate formation kinetic models is found in Vysniauskas & Bishnoi (1983, 1985). They proposed a semi-empirical model based on the Arrhenius equation to describe the formation of critical methane and ethane hydrate clusters and their subsequent growth in a semi-batch stirred tank reactor. Following their work, similar studies have been conducted to define the growth kinetics of gas mixtures and single gas hydrates. These investigations either adopt new apparatus, scale up the analysis, introduce new parameters for consideration, or provide alternative approaches (Al-Otaibi et al., 2010; Chun & Lee, 1996; Clarke & Bishnoi, 2005; Englezos et al., 1987; Lekvam & Ruoff, 1993). Pioneering studies on mass transfer mechanisms were carried out by Herri et al. (1999) and Skovborg & Rasmussen (1994). Englezos' model of intrinsic kinetics was modified by Hashemi et al. (2007) by replacing the driving force of fugacity difference with concentration difference. It argued that the overall kinetic rate constant is difficult to obtain, while the growth process is more controlled by mass diffusion. This forms the basis of subsequent mass transfer-limited
models, as seen in Turner et al. (2009) and Mohebbi et al. (2012). The surrounding medium absorbs the heat released from an exothermic reaction. The exothermic nature of hydrate formation activated the study to include this effect on the overall growth of hydrate molecules. (Uchida et al. 1999) carried out pioneering research on this effect, concluding that the linear growth rate of the hydrate film depends more on the degree of subcooling than on pressure. Subsequent studies focused on concerns about the impact of counter-current convection, convective heat transfer, backward heat transfer, and conductive heat transfer (Freer et al., 2001; Mochizuki & Mori, 2006; Mori, 2001).
In the growing field of hydrate kinetics, the sparse available growth models have not been inclusive of the three known mechanisms that control the growth of gas hydrates, creating a need in the industry for a more comprehensive model. The advantages of a comprehensive study that yields a complete model are that it enables the mapping of a more accurate kinetic growth profile representation of gas hydrates. The application of previous models is limited. Some models present kinetic growth via intrinsic kinetics only, other models may show growth as a function of heat evolution and mass transfer, only overlooking the impact of intrinsic kinetics. As the transportation system conditions can change, they are most likely to encounter conditions that may favour the dominance of one controlling gas hydrate growth mechanism over another. Without a comprehensive model, it will be challenging to switch from one model to another when a particular dominating mechanism is suspected. A complete model makes is easier to define the growth prfile of hydrate in any gas transportation system. As an example, For systems such as a turbulent multiphase pipeline flow, intrinsic kinetics may play only a minor yet significant role in the formation of hydrate, while mass and heat transfer would to a large extent determine the rate of hydrate formation (Sloan & Koh, 2007; Meindinyo, 2017; Ke et al., 2019); in high-production gas lines (where influence of mass and heat transfer become negligible), shutdown operations (with a dominating mass and heat transfer effect), and insulated pipelines with increased heat transfer from hydrate growth. Recent hydrate growth studies by Ke et al. (2019), Yin et al. (2018), and Filarsky et al. (2023) have also highlighted the need for a comprehensive model that will incorporate the controlling mechanisms of heat transfer, mass transfer and intrinsic kinetics. This provides a flexible tool that can be employed in quantifying hydrate growth in a variety of pipeline system conditions.
-
METHODOLOGY
-
Advection-Diffusion Transport Equation
The advection-diffusion equation is a combination of the diffusion and advection equations, describing a physical phenomenon in which particles, energy, or other physical quantities are transferred within a physical system due to diffusion and advection. Advection is the transport of a substance or quantity by bulk motion. At the same time, diffusion relates to the movement of a substance from an area of high concentration to an area of low concentration. The advective diffusive equation is fundamentally derived from the continuity equation, which states that in a differential control volume, the rate of change of a scalar quantity is given by the net flow and diffusion into and out of that part of the system, considering any generation or consumption within the control volume. The ADE has been applied in the studies of wax deposition in subsea pipelines (Lee, 2008; Huang, 2011) and for methane hydrate studies in porous media (Rempel & Buffett, 1998; Hesse, 2003; Haeckel et al., 2004; Peszynska et al., 2015; VanderBeek & Rempel, 2018). These studies show that the advective-diffusive transport equation may be modified to account for the hydrate growth rate within a control volume, taking into account the influence of quantities such as heat, mass, and intrinsic kinetics.
The proposed model improves on the model of Freer et al. (2001) into the ADE. The influence of mass concentration gradients forms the base of the model while heat production and dissipation, and the intrinsic reaction property of the hydrate growth process is modelled into the reactive component of the equation. The basic form of the equation is given in Equation 1 as (Miyaoka et al., 2017):
+ ((, ) + (, )) = (, , , ) (1)
Where, for mass flux within the control volume, denotes the concentration of hydrate particles, we calculate the Laplacian () operator in spatial variables and ; (, ) is the diffusion coefficient, (, ) = ((, ), (, )) is the velocity of the flowing fluid; and (, , , ) is an external concentration source, which includes reaction, dynamic terms, depending on or not on . The quantity may also represent heat flux in the heat transfer ADE.
Figure 1: Schematic illustration of film growth at the gas-liquid interface. A moving boundary model adapted from Freer et al. (2001)
The velocity of the fluid in the x-direction (see Figure 1) results from density changes at the moving boundary
-
The proposed expression for hydrate growth is given by:
Where:
+
=
2
2
( ) (2)
1 = 1 + 1
(2b)
=
(2c)
= (
)
(2d)
-
IC: (, 0) = , 0 0
-
BC 1: (0, ) =
-
BC 2: (, ) = > 0, = ()
Where is the intrinsic kinetic constant, and is the heat transfer coefficient.
-
-
Analytical and numerical models
To translate the kinetic growth of methane gas hydrates resulting from the three driving mechanisms into solution terms, both analytical and numerical modelling approaches are attempted. The analytical model makes use of the separation of variables technique in solving the proposed kinetic model. The method of separation of variables is a powerful technique used to solve partial differential equations (PDEs). This method involves breaking down a complex PDE into simpler ordinary differential equations (ODEs) by assuming that the solution can be expressed as a product of functions, each depending on only one variable. The numerical model utilises Taylors series second-order expansion to discretise the proposed model, making it suitable for numerical evaluation and implementation. The Crank-Nicolson (CN) finite difference method is then used to solve the discretised model. John Crank and Phyllis Nicolson developed the CN method as a numerical solution of a PDE which arises from heat conduction problems (Crank & Nicolson, 1996). It was introduced to curb instability and to increase the efficiency and accuracy of both implicit and explicit methods. The choice of this method is based on the fact that its accuracy is of second order in space-time discretisation, as well as unconditional stability in time (Umeorah & Mashele, 2019). In discretising, two significant challenges are faced: stability and accuracy. Stability is improved using a range of implicit methods while accuracy requires the best form of spatial discretisation.
In discretising our model equation, implicit schemes are used for their stability and Taylor expansion for their accuracy in the spatial coordinate. Since the Crank-Nicolson method combines both the explicit Euler FTCS and the implicit Euler BTCS, its accuracy is higher than the two aforementioned finite difference methods.
Figure 2: The Mesh for Crank Nicholson Scheme
It should be noted that the Crank-Nicolson methods are numerically unconditionally stable (Umeorah & Mashele, 2019). A typical mesh for the scheme is shown in Figure 2.
-
-
Model solutions
Mathematical modelling plays a pivotal role in understanding physical, biological, and engineering phenomena. Many such models involve partial differential equations (PDEs), but analytical solutions are not always feasible. When possible, analytical methods provide exact solutions, while numerical techniques offer approximations when closed-form solutions are unavailable.
-
Analytical solution
The model may be expressed as:
2
= 2 + (3)
Where is regarded as a constant, in this case, the sum of intrinsic kinetics and heat transfer functions. Solving the equation above (Equation 3) requires first solving the partial differential equation and extending the solution to include . The advection-diffusion equation to solve becomes:
2
= 2 (4)
We first reduce the ADE in Equation 4 into a diffusion equation to solve this analytically. This is achieved by defining a function y such that:
=
= 1,
=
(, ) = (, )
Defining with respect to the new variable by applying the chain rule:
=
2 2
2 = 2
Since =
= +
=
Substituting the derivatives into the differential equation, we get:
2
=
2
2
= 2 (5)
The convection terms on both sides cancel out, and the equation is reduced to the diffusivity equation (Equation 5). We will need to transform the initial condition as well, but and are equal when = 0
(, 0) = (, 0)
Now we try to solve the diffusivity equation. However, the problem presented in the model has different specified boundary conditions (inhomogeneous boundary conditions):
(0, ) = 1, (, ) = 2, > 0.
The concentration (, ) satisfies the diffusion equation with diffusivity :
2
= 2
=
Approaching this problem directly with the separation of variables could prove tricky and may not yield results. Therefore, we will consider alternative methods for a more effective solution. For this solution, let us use as for simplification, then applying the inhomogeneous boundary condition at = 0 directly to the ansatz (, ) = ()() results in:
(0, ) = (0)() = 1;
So that,
(0) = 1/().
However, our separation of variables ansatz assumes () to be independent of . We can therefore say that inhomogeneous boundary conditions are not separable. The proper way to solve a problem with inhomogeneous boundary conditions like we have is to transform it into another problem with homogeneous boundary conditions. As , we assume that a stationary concentration distribution () will attain, independent of . Since () must satisfy the diffusion equation, we have:
"() = 0, 0
With a general solution:
() = + .
Since () must satisfy the same boundary conditions of (, ), we have (0) = 1 and () = 2, and we determine = 1 and
= (2)/. We express (, ) as the sum of the known asymptotic stationary concentration distribution () and an unknown transient concentration distribution (, ):
(, ) = () + (, ).
Substituting into the diffusion equation, we obtain
2
(() + (, )) = 2 (() + (, ))
Or
= ,
Since = 0 and = 0. The boundary conditions satisfied by are:
(0, ) = (0, ) (0) = 0,
(, ) = (, ) () = 0,
So that is observed to satisfy homogeneous boundary conditions. If the initial conditions are given by (, 0) = (), then the initial conditions for are:
(, 0) = (0, ) (0) = 0
= () ().
The resulting equations may then be solved for (, ) using the technique for homogeneous boundary conditions, and (, ) subsequently determined. Then, introducing the homogeneous boundary condition technique, we write the 1 diffusion problem in a domain of length , and solve the diffusion equation for the concentration (, ),
2
= 2 , 0 , > 0.
Both initial and boundary conditions are required for a unique solution. That is, we assume the initial concentration distribution in the pipe section of interest is given by:
(, 0) = (), 0 .
Furthermore, we assume that boundary conditions are given at the ends of the pipes. When the concentration value is specified at the boundaries, the boundary conditions are called Dirichlet boundary conditions. We assume here homogeneous Dirichlet boundary conditions, that is zero concentration at both boundaries, which could occur if the ends of the medium open up to into large reservoirs of clear solutions.
(0, ) = 0 (, ) = 0, > 0.
If () is identically zero, then the trivial solution (, ) = 0 satisfies the differential equation and the initial and boundary conditions and is therefore the unique solution of the problem. In what follows, we will assume that () is not identically zero, so we need to find a solution different from the trivial solution.
The method we utilize is referred to as the separation of variables, a powerful technique often used to solve partial differential equations. In this approach, we assume that the solution can be expressed as the product of two distinct functions: one function that varies only with respect to position and another function that varies solely concerning time. This separation allows us to isolate the variables, making it possible to analyse each component individually and ultimately solve the equation more effectively. By applying this technique, we can transform a complex problem into simpler, more manageable parts, facilitating a clearer path to finding the solution.
(, ) = ()()
()
=
2()
2
Because is only a function of position and is only a function of time,
2
= 2
1
=
Where is the second spatial derivative of and is the time derivative of .We then assume that the above equation is equal to a constant 2, such that:
= 1 = 2
Now we have succeeded in reducing our initial problem into two ODEs
= 2
Therefore,
1 = 2
+ 2 = 0 1
+ 2 = 0 2
First, we solve 2 using a set of boundary equations to get and then substitute into 1. Because is second derivative, we need two boundary conditions to solve , and if first derivative needs one condition.
+ 2 = 0
Has the general solution:
() = +
Applying the boundary conditions:
(0, ) = 0
Therefore, (0)() = 0. For this to hold, (0) = 0
(, ) = 0
Therefore, ()() = 0. For this to hold, () = 0.
Putting (0) = 0 into the general equation, (0) = 0 = . Therefore, = 0 (0 = 0, 0 = 1).
() = 0. Therefore, = 0 (Since is already 0). is a constant, so = 0. will be equal to zero when =
for = 1,2,3,
Therefore,
And
Is the solution to 2.
= /
() =
For 1, we substitute 2 = 22 into the equation.
2
Let 22 = 2 and
=
+
22
2 = 0
2
+ 2 = 0
Transforming the expression and dividing both sides by :
= 2
1 = 2 0
Integrating both sides and evaluating the integrals,
1 = 2
ln = 2 +
Converting the logarithm into exponential form
() =
2 +
=
2
×
() =
2
But (, ) = ()(). Now we put together the solutions to the two ODEs obtained:
(, ) =
2 ×
sin
Where × = ,
2
(, ) = sin
And,
2
(, ) = sin
=1
The final solution step is to satisfy the initial conditions to determine what the constant represents. Therefore, at = 0, we have
(, ) = ()
() = sin
=1
2
= () sin
0
(6)
Given () = 1, the integral of 1 sin / over [0, L] is:
=
2
sin
=
21 [
cos
]
= 21 (1 (1))
1
0
0
(, ) = 2 sin
But,
=1
(, ) = () + (, ).
2 2
(, ) = + + sin (7)
1
=1
(, ) is a solution to the diffusivity equation (Equation 5). Rewriting this solution in terms of which was used as , it becomes:
(, ) =
+ 2 +
sin 2
1
=1
Substituting = into this solution above gives the complete analytical solution to the advection- diffusion equation (Equation 4).
(, ) =
+ 2( ) +
sin ( ) 2 (8)
1
=1
And the solution to our comprehensive model (Equation 3) becomes:
2( ) ( ) 2
(, ) = + + sin + (9)
1
=1
The modelled solution formula can be interpreted to include: a constant concentration, 1; a linear advection component,
2(); a transient decaying term (sine series); and an offset, that absorbs the intrinsic kinetics and heat constraints.
-
Numerical Solution
The equation:
2
( )
= 2 +
Let ()
be equal to constant as in Equation 3. Using the C-N scheme (which is implicit and
second-order accurate in time), we approximate:
+1
=
And take an average of spatial terms at time levels and + 1. This gives:
2 1
+1
+1
+1
2 = 2()2 {(+1
2
+ 1 ) + (+1 2
+ 1)}
1 +1 +1
1
= ( +1 1 + +1 1) =
(+1 +1 + )
2
2
2
4
+1
1
+1
1
Now, substituting these approximations into the original PDE, the complete discretised equation looks like:
1
+1
= [
{(+1 2+1 + +1) + ( 2 + )}]
2()2
+1
1
+1
1
1
[
(+1 +1 + )]
Let =
22
=
4
4
+1
1
+1
1
+1 = [(+1 2+1 + +1) + ( 2 + )] [(+1 +1 + )]
+1
1
+1
1
+1
1
+1
1
Resolving the unknown j+1 to the left and known j to the right, the equation becomes
+1 (+1 2+1 + +1) + (+1 +1) = + ( 2 + ) ( ) +
+1
1
+1
1
+1
1
+1
1
Further simplification of the discretised model gives
( )+1 + (1 + 2)+1 ( )+1 = ( + ) + (1 2) + ( ) +
1
+1
1
+1
The equation makes sense for space indices = 1, , 2 but it does not make sense for indices on the boundaries: = 0 and = 1.
= 0: ( )+1 + (1 + 2)+1 ( )+1 = ( + ) + (1 2) + ( ) +
1 0 1 1 0 1
= 1: ( )+1 + (1 + 2)+1 ( )+1
2 1
= ( + ) + (1 2) + ( ) +
2 1
But and lie outside the grid beyond the boundary. Therefore, considering Neumann boundary
1
condition,
= 0: = +1 = +1
0 1 0 1
= 1: = +1 = +1
1
= [
0
]
1
1
Using this notation, the above approximation for a fixed point in time, = , can be written compactly as a linear system:
(1 + ) +
0 0 0
0 0 0
+1
0
0
+1
1 + 2 + 0 0 0 0 0 0 1
0
1 + 2
+
0 0 0 0
0 +1
0 0
0 0 0
0 0 0 2
0 0 0 1 + 2 + +1
[ 0 0 00 0 0
0 1 + + ]
2
(1 + )
+ 1 2
0 0 0
0 0
0 0 0
0 0 0
+1
1
0
0
[ ]0
1
= 0 +
0 0
1 2
0 0 0 0
0 0
0
2
0
0 0 0
0 0 0 +
1 2
[ 0 0 00 0 0 0
+
1 ]
2
[ ]+
[]
The matrix above is in the form:
+1 = +
1
Where represents the constant on the RHS of the equation. Therefore, the numerical solution becomes:
+1 = 1( + )
Stability Analysis
We employ the Von Neumann stability analysis for our Crank-Nicolson scheme, The numerical scheme is given as:
( )+1 + (1 + 2)+1 ( )+1 = ( + ) + (1 2) + ( )
1
+1
1
+1
We assume a Fourier mode solution of form:
=
Where is the amplification factor, is the wave number, is the spatial step in this case. Substituting this into the numerical scheme:
= (+1) =
+1
= (1) =
1
Similarly for the time step + 1:
+1 = +1+1
+1
+1 = +1+1
1
Rewriting the equation in terms of Fourier modes:
( )+1 + (1 + 2)+1 ( )+1
= ( + ) + (1 2) + ( )
Now dividing by and setting +1 = , will give:
( ) + (1 + 2) ( ) = ( + ) + (1 2) + ( )
To solve for , we define cos() and sin():
+ = 2 cos()
= 2 sin()
Using these definitions, we rewrite:
( )() + (1 + 2) ( )() = ( + ) + (1 2) + ( )
Grouping:
[(1 + 2) ( ) ( + )] = (1 2) + ( ) + ( + )(1 2) + ( ) + ( + )
= (1 + 2) ( ) ( + )
(1 2) + 2 cos() ( )
= (1 + 2) 2 cos() ( + )
For stability, we require: || 1
For unconditional stability, || 1 for all values of . This means:
|(1 2) + 2( ) cos()| |(1 + 2) 2( + ) cos()|
This must hold for all values, in the worst case, when cos() ± 1: For cos() = 1,
This is marginally stable. For cos() = 1,
(1 2) + 2( )
= (1 + 2) 2( + ) 1 + 2 2
= = 1
1 + 2 2
For unconditional stability,
(1 2) 2( )
= (1 + 2) + 2( + )
1 4 + 2
=
1 + 4 + 2
For stability:
1 4 + 2
|
1 + 4 + 2
| 1
1
1 4 + 2
1
1 + 4 + 2
If we take the first inequality and multiply by the denominator, we get
(1 + 4 + 2) 1 4 + 2
Solving for , 1.
2
If we take the second inequality:
1 4 + 2
1
1 + 4 + 2
Multiplying by the numerator and solving for , 0.
For unconditional stability, the minimum required values are 1 and 0. This means that must be
2
non-negative and must be at least 1.
2
-
-
SIMULATION
-
Analytical Simulation result
We present analytical results for the model solution derived in the previous chapter. In generating these results, we use the following values for parameters: constants 1 and 2 derived from gas equations.
-
Hydrate concentration propagation over time
The analysis of our analytical model results shows that the concentration plot of methane gas hydrates formed by the combined influence of all three mechanisms shifts to the right and exhibits a higher increase in time difference. This is illustrated in Figure 3. This shift towards the right tends to produce lower values for concentration at specific spatial values.
he plotted graph in Figure 3 shows how the hydrate concentration, (, ) evolves over space [0,1] for various time points = 0, 500, 1000, 5000. The following trend is observed at the different times: At = 0, The solution includes strong spatial oscillations due to the full presence of the Fourier sine terms. The concentration is at its most non-uniform state, reflecting the initial sharp variation (like an impulse or discontinuity in the initial condition). As increases (500, 1000, 5000): The oscillations smooth out, and the solution gradually transitions to a more linear profile. This is caused by the exponential decay of higher-order modes
(2 ), effectively damping the sine terms. The concentration approaches a steady-state linear profile: (, ) +
1
2( ) + which represents diffusion-dominated transport. Long-term behaviour (e.g., = 5000): The system is nearly fully
diffused, and the impact of the initial Fourier series has almost vanished. The profile becomes smooth and slowly shifts due to advection (the term).
Figure 3: Hydrate concentration propagation over time
-
Influence of Advection
The advection velocity shifts the spatial position in time. This causes the entire concentration profile to shift over time, moving in the direction of the flow. This is illustrated in Figure 4.
Figure 4: Effect of advection of the hydrate growth profile
As increases, decreases (for > 0) causing the rightward shift over time is experienced. This shift represents the movement of the hydrate concentration due to advection represented by . The speed of this shift is controlled by the magnitude of . As can be seen, advection does not directly change the shape of the solution but rather repositions it over time.
-
Influence of diffusion
Diffusion causes the smoothing or flattening of the initial distribution. It acts to harmonise the concentration profile by spreading mass from high to low regions.
Figure 5: Effect of diffusion on hydrate growth profile
The plot (Figure 5) illustrates the influence of the diffusion coefficient on the hydrate concentration profile (, ) at a fixed time
= 1000. At low diffusion between the orders 1018 to 1010, the concentration profile tends to retain the sharp variations and oscillatory features from the initial condition. Very little smoothing occurs which translates to the system dominated by advection. The profile evolves slowly and remains close to the original configuration for small . The profile begins to smooth out between orders of = 108 to 105, peak values decrease and valleys fill in to indicate the decay of high-frequency modes. For higher diffusion (at = 103), the profile becomes nearly linear and smooth across . Strong damping of all but the most slowly varying
modes results in near steady-state behaviour. The shape resembles the background linear profile 1
2
+ ( )
( ) + , with little
contribution from the sine terms. The influence of the initial condition is almost lost and the solution is diffusion dominated. The interaction of advection and diffusion creates the dynamic spreading and shifting behaviour seen in the plot hydrate concentration over time.
-
Influence of the initial condition
The initial condition represents a key parameter that can shape the hydrate concentration profile in the derived model. Initial condition (via Fourier terms) strongly affects the early-time behaviour but quickly decays for large . The impact is illustrated in Figure 6.
Figure 6: Effect of initial condition on the hydrate profile
As 1 increases, the entire concentration profile shifts upward. This is expected as 1 contributes both a constant term and scales the Fourier sine coefficients , which affects the transient part of the solution. The overall shape of the profile remains similar
across different 1 values showing that the initial concentration magnitude affects the amplitude, but not the structure of the spatial distribution at time . At time = 1000 the effect of the initial sine components is reduced due to diffusion, however, because the sine terms scale with 1, larger 1 still leads to visibly higher amplitude undulations. Therefore, higher 1 results in a higher mean concentration and a stronger (though still decaying) residual transient effect from the initial condition.
-
Influence of the M parameter
In the solution equation (Equation 9), the parameter acts as a constant additive offset. Increasing uniformly shifts the entire concentration curve upwards. The shape, slope, and spatial variation of the profile remain unchanged.
Figure 7: Effect of intrinsic kinetics and heat transfer on the hydrate growth profile
There is no impact on the plot's dynamics as seen in Figure 7. Since is time-independent (in the solution equation) and not involved in the spatial derivatives or decay terms, it does not affect how the profile evolves, only its vertical position. That means diffusion, advection, and initial condition dynamics remain unaffected. denotes the baseline concentration, which encapsulates the effects of intrinsic kinetics associated with the initial growth reaction, as well as the limitations imposed by heat transfer. This concentration serves as a foundational measure that integrates the rate of reaction occurring within the system and the efficiency of heat dissipation or retention, both of which are crucial in determining the overall dynamics of the process. Understanding is essential for analysing how these factors interact to influence the performance and behaviour of the system.
-
-
Numerical simulation results
-
Hydrate concentration propagation over time
Figure 8: Numerical solution profile transition with time
Figure 8 shows the spatial concentration profile at times between zero and 5000 seconds. Initially, the entire domain is represented by the initial concentration and the base concentration (intrinsic kinetics and heat constraints), which shows a flat line. As time progresses, a gentle gradient forms across the domain; this gradient tends to reflect a slow diffusion from the higher concentration at the right boundary towards the left. Advection is very weak due to the small so the shape is primarily governed by diffusion. At times 500s and 1000s the profile starts bending upward slightly on the right but remains nearly flat. At 5000s, there is a noticeable but still slow change; the system is yet to reach a steady state due to very small diffusion and very slow advection. This very slow transport in the physical system is a representation of molecular diffusion in gas hydrate buildup.
-
Influence of diffusion on the numerical solution profile
This sensitivity analysis is performed by varying the diffusion coefficients in the solution analysis and observing the trend as they transition in time. A plot of the analysis is provided in Figure 9.
In Figure 9, the four sub-plots correspond to different diffusion coefficients between = 1016 and 104 showing concentration profiles progressing from = 0 to 5000. At very low diffusion, where = 1016, the profile changes very slowly and almost flat at 5000s showing minimal diffusion effects. At the times shown, the higher boundary has barely influenced the interior. As the diffusion coefficient increases to 1010, there is a slight smoothing of the profile over time and gradients begin to form although the centre is still mostly flat at early times but some visible diffusion by 5000. On further increase of to 106, a clear gradient from left boundary to right boundary forms faster and the effect of diffusion is observed earlier. At 5000s, profile is more uniformlytransitioning between boundaries. On increasing the value of to 104, we see that diffusion quickly begins to dominate.
Figure 9: Sensitivity of Diffusivity in the numerical solution profile
There is a quick transition to a smooth curve between boundary values and a nearly linear gradient between left and right ends early in the simulation. A higher diffusion coefficient spreads concentration differences more rapidly and moves in the direction
of flow for positive advection. However, when the coefficient is low, the profile will largely depend on other controlling mechanisms and take a longer time for the diffusion effect to be observed.
-
Influence of on the Numerical Solution Profile
This sensitivity analysis is performed by varying the values in the solution analysis and observing the trend as they transition in time. A plot of the analysis is provided in Figure 10.
Figure 10: Effect of on the numerical solution profile
Figure 10 shows the concentration profiles with and without the offset that represents the intrinsic kinetics and heat transfer constraints. It clearly illustrates that the shape of the profiles remains unchanged as does not impact the physics of the simulation. Also, the magnitude is uniformly increased by , indicating a vertical shift upwards.
-
The initial condition on the Numerical Solution Profile
This sensitivity analysis is performed by varying the initial conditions in the solution analysis and observing the trend as they transition in time. A plot of the analysis is provided in the Figure 11. The plot shows the numerical solution profiles at values lower and higher than the baseline value. For each value, the concentration profile (, ) is shown at different times, illustrating how increasing 1 raises the entire profile and enhances early transport towards the right boundary. The system is linearly sensitive to changes in the initial condition. This linearity reflects the boundary-dominated nature of the solution in the low-diffusion regime. Because the diffusion is so small, information propagates very slowly, and the initial boundary conditions that include the preexisting controlling mechanism, dominate for a while.
Figure 11: Effect of initial condition on the numerical solution profile
-
-
-
Discussion of solution results
-
In solving the one-dimensional advectiondiffusion equation for hydrate growth, with inhomogeneous boundary condition, a consistent discrepancy was observed between the analytical and numerical solutions. This disagreement may be attributed to a number of factors intrinsic to numerical simulation. Finite-difference schemes of second-order accuracy, although widely adopted, inevitably introduce truncation errors that can accumulate during integration. While coarse spatial discretization is a common source of error, the present work employed sufficiently fine grid spacing, making this effect unlikely to dominate. The CrankNicolson method, used for time advancement, is theoretically unconditionally stable; stability of the implementation was verified for the parameter space considered, thus reducing the likelihood of instability-related deviations.
Residual differences may arise from imperfect replication of the idealized boundary and initial conditions assumed in the analytical derivation. Discrete representations cannot perfectly reproduce continuous profiles, and even small mismatches can propagate over time. In advection-dominated regimes, numerical dispersion and phase errors may further distort wave propagation, leading to shifts or attenuation relative to the analytical form. Such artefacts can also be introduced by averaging operations implicit in the discretization, which may smooth gradients in a manner inconsistent with the analytical formulation.
Additional discrepancies may stem from modelling assumptions. Analytical solutions typically rest on simplified conditionssuch as spatial homogeneity, constant transport coefficients, and negligible boundary-layer effectsthat are not always replicated in numerical implementation. Moreover, computational errors, including incorrect indexing in discretization stencils, misapplication of boundary conditions, and solver inaccuracies, can bias results. Rounding errors from floating-point arithmetic, though generally minor, may accumulate in long-duration simulations.
Addressing these discrepancies requires systematic verification and validation and is beyond the scope of this work. This includes grid and time-step refinement studies to confirm convergence rates, careful alignment of numerical boundary and initial conditions with the analytical case, and benchmarking against manufactured or canonical solutions. Through such measures, numerical models can be refined to reduce deviations, thereby enhancing both the reliability and predictive capacity of hydrate growth simulations.
Conclusion and recommendation
A comprehensive novel kinetic model for methane hydrate growth in gas flowlines was developed. In this model, the controlling mechanisms responsible for gas hydrate growth are represented. This is the first work among published literature that takes into account all three mechanisms. Two solutions methods are then presented for the proposed comprehensive model for methane gas hydrates. Analytical solution solved using the separable variable technique, and the Numerical solution using the Crank- Nicolson method. Results indicate a similar trend exhibited by solution framework parameters. However, within the numerical and analytical solution results, there exist discrepancies in the magnitude of the concentration profile. This inconsistency may be a result of a number of factors, from solution implementation to assumptions application.
Understanding the sources of discrepancies between the numerical and analytical solutions of our proposed comprehensive hydrate growth equation is crucial for further reasearch. In addressing these factors, researchers can refine the numerical model in the future to achieve better agreement with the analytical results, leading to more reliable and accurate simulations. Further
validation of numerical methods against well-understood systems can also help enhance their credibility and applicability. In the quest for better prediction and management of natural gas hydrates, future research can dwell on upgrading hydrate simulation software like OLGA and HWPVT or developing algorithms that can integrate or couple the thermodynamic models together with the kinetic growth model. This way, the occurrence of hydrate can be predicted and quantified by a singular software.
REFERENCES
-
Al-Otaibi, F., Clarke, M., Maini, B., & Bishnoi, P. R. (2010). Formation kinetics of structure i clathrates of methane and ethane using an in situ particle size analyzer. Energy and Fuels, 24(9), 50125022. https://doi.org/10.1021/EF100560F
-
Arjun, A., & Bolhuis, P. G. (2020). Rate Prediction for Homogeneous Nucleation of Methane Hydrate at Moderate Supersaturation Using Transition Interface Sampling. Journal of Physical Chemistry B, 124(37), 80998109. https://doi.org/10.1021/ACS.JPCB.0C04582/SUPPL_FILE/JP0C04582_SI_003.MP4
-
Chun, M. K., & Lee, H. (1996). Kinetics of formation of carbon dioxide clathrate hydrates. Korean Journal of Chemical Engineering, 13(6), 620626. https://doi.org/10.1007/BF02706029/METRICS
-
Clarke, M. A., & Bishnoi, P. R. (2005). Determination of the intrinsic kinetics of CO2 gas hydrate formation using in situ particle size analysis. Chemical Engineering Science, 60(3), 695709. https://doi.org/10.1016/J.CES.2004.08.040
-
Crank, J., & Nicolson, P. (1996). A practical method for numerical evaluation of solutions of artial differential equations of the heat-conduction type. Advances in Computational Mathematics, 6(1), 207-226.
-
Creek, J. L., Subramanian, S., & Estanga, D. A. (2011). New Method for Managing Hydrates in Deepwater Tiebacks. 25. https://doi.org/10.4043/22017-MS
-
Englezos, P., Kalogerakis, N., Dholabhai, P. D., & Bishnoi, P. R. (1987). Kinetics of formation of methane and ethane gas hydrates. Chemical Engineering Science, 42(11), 26472658. https://doi.org/10.1016/0009-2509(87)87015-X
-
Filarsky, F., Hagelstein, M., & Schultz, H. J. (2023). Influence of different stirring setups on mass transport, gas hydrate formation, and scale transfer concepts for technical gas hydrate applications. Applied Research, 2(1), e202200050.
-
Freer, E. M., Sami Selim, M., & Dendy Sloan, E. (2001). Methane hydrate film growth kinetics. Fluid Phase Equilibria, 185(12), 6575. https://doi.org/10.1016/S0378-3812(01)00457-5
-
Haeckel, M., Suess, E., Wallmann, K., & Rickert, D. (2004). Rising methane gas bubbles form massive hydrate layers at the seafloor. Geochimica et Cosmochimica Acta, 68(21), 4335-4345.
-
Hashemi, S., Macchi, A., & Servio, P. (2007). Gas hydrate growth model in a semibatch stirred tank reactor. Industrial and Engineering Chemistry Research, 46(18), 59075912. https://doi.org/10.1021/IE061048
-
Herri, J. M., Pic, J. S., Gruy, F., & Cournil, M. (1999). Methane hydrate crystallization mechanism from in-situ particle sizing. AIChE Journal, 45(3), 590 602. https://doi.org/10.1002/AIC.690450316
-
Hesse, R. (2003). Pore water anomalies of submarine gas-hydrate zones as tool to assess hydrate abundance and distribution in the subsurface: What have we learned in the past decade?. Earth-Science Reviews, 61(1-2), 149-179.
-
Huang, Z. (2011). Application of the fundamentals of heat and mass transfer to the investigation of wax deposition in subsea pipelines (Doctoral dissertation, University of Michigan).
-
Kashchiev, D., & Firoozabadi, A. (2002a). Driving force for crystallization of gas hydrates. Journal of Crystal Growth, 241(12), 220230. https://doi.org/10.1016/S0022-0248(02)01134-X
-
Kashchiev, D., & Firoozabadi, A. (2002b). Nucleation of gas hydrates. Journal of Crystal Growth, 243(34), 476489. https://doi.org/10.1016/S0022- 0248(02)01576-2
-
Ke, W., Svartaas, T. M., & Chen, D. (2019). A review of gas hydrate nucleation theories and growth models. Journal of Natural Gas Science and Engineering, 61, 169196. https://doi.org/10.1016/J.JNGSE.2018.10.021
-
Kelkar, S. K., Selim, M. S., & Sloan, E. D. (1998). Hydrate dissociation rates in pipelines. Fluid Phase Equilibria, 150151(151), 371382. https://doi.org/10.1016/S0378-3812(98)00337-9
-
Lee, H. S. (2008). Computational and Rheological Study of Wax Deposition and Gelation in Subsea Pipelines (Doctoral dissertation).
-
Lekvam, K., & Ruoff, P. (1993). A Reaction Kinetic Mechanism for Methane Hydrate Formation in Liquid Water. Journal of the American Chemical Society, 115(19), 85658569. https://doi.org/10.1021/JA00072A007
-
Max, M. D., Johnson, A. H., & Dillon, W. P. (2006). Economic Geology of Natural Gas Hydrate. 9. https://doi.org/10.1007/1-4020-3972-7
-
Meindinyo, R. E. T. (2017). Gas Hydrate Growth Kinetics: Experimental Study Related to Effects of Heat Transfer.
-
Miyaoka, T. Y., Meyer, J. F. D. C. A., & Souza, J. M. R. (2017). A general boundary condition with linear flux for advection-diffusion models. TEMA (São Carlos), 18(2), 253-272.
-
Mochizuki, T., & Mori, Y. H. (2006). Clathrate-hydrate film growth along water/hydrate-former phase boundaries-numerical heat-transfer study. Journal of Crystal Growth, 290(2), 642652. https://doi.org/10.1016/j.jcrysgro.2006.01.036
-
Mohebbi, V., Naderifar, A., Behbahani, R. M., & Moshfeghian, M. (2012). Investigation of kinetics of methane hydrate formation during isobaric and isochoric processes in an agitated reactor. Chemical Engineering Science, 76, 5865. https://doi.org/10.1016/j.ces.2012.04.016
-
Mori, Y. H. (2001). Estimating the thickness of hydrate films from their lateral growth rates: Application of a simplified heat transfer model. Journal of Crystal Growth, 223(12), 206212. https://doi.org/10.1016/S0022-0248(01)00614-5
-
Peszynska, M., Showalter, R. E., & Webster, J. T. (2015). Advection of methane in the hydrate zone: model, analysis and examples. Mathematical Methods in the Applied Sciences, 38(18), 4613-4629.
-
Rempel, A. W., & Buffett, B. A. (1998). Mathematical models of gas hydrate accumulation. Geological Society, London, Special Publications, 137(1), 63-74.
-
Skovborg, P., & Rasmussen, P. (1994). A mass transport limited model for the growth of methane and ethane gas hydrates. Chemical Engineering Science, 49(8), 11311143. https://doi.org/10.1016/0009-2509(94)85085-2
-
Sloan Jr., E. D., & Koh, C. A. (2007). Clathrate Hydrates of Natural Gases. Clathrate Hydrates of Natural Gases. https://doi.org/10.1201/9781420008494
-
Stoporev, A. S., Semenov, A. P., Medvedev, V. I., Kidyarov, B. I., Manakov, A. Y., & Vinokurov, V. A. (2018). Nucleation of gas hydrates in multiphase systems with several types of interfaces. Journal of Thermal Analysis and Calorimetry, 134(1), 783795. https://doi.org/10.1007/S10973-018-7352-2/FIGURES/8
-
Turner, D. J., Miller, K. T., & Dendy Sloan, E. (2009). Methane hydrate formation and an inward growing shell model in water-in-oil dispersions. Chemical Engineering Science, 64(18), 39964004. https://doi.org/10.1016/J.CES.2009.05.051
-
Uchida, T., Ebinuma, T., Kawabata, J., & Narita, H. (1999). Microscopic observations of formation processes of clathrate-hydrate films at an interface between water and carbon dioxide. Journal of Crystal Growth, 204(3), 348356. https://doi.org/10.1016/S0022-0248(99)00178-5
-
Umeorah, N., & Mashele, P. (2019). A Crank-Nicolson finite difference approach on the numerical estimation of rebate barrier option prices. Cogent Economics & Finance.
-
VanderBeek, B. P., & Rempel, A. W. (2018). On the importance of advective versus diffusive transport in controlling the distribution of methane hydrate in heterogeneous marine sediments. Journal of Geophysical Research: Solid Earth, 123(7), 5394-5411.
-
Vysniauskas, A., & Bishnoi, P. R. (1983). A kinetic study of methane hydrate formation. Chemical Engineering Science, 38(7), 10611072. https://doi.org/10.1016/0009-2509(83)80027-X
-
Vysniauskas, A., & Bishnoi, P. R. (1985). Kinetics of ethane hydrate formation. Chemical Engineering Science, 40(2), 299303. https://doi.org/10.1016/0009- 2509(85)80070-1
-
Yin, Z., Chong, Z. R., Tan, H. K., & Linga, P. (2016). Review of gas hydrate dissociation kinetic models for energy recovery. Journal of Natural Gas Science and Engineering, 35, 1362-1387.
-
Yin, Z., Khurana, M., Tan, H. K., & Linga, P. (2018). A review of gas hydrate growth kinetic models. Chemical Engineering Journal, 32, 929. https://doi.org/10.1016/J.CEJ.2018.01.120
