Mathematical Modeling of Wax Deposition in Field-Scale Crude Oil Pipeline Systems

. The formation of solid wax crystals, which interlock and form a gel-like layer on the inner wall of the crude oil pipeline, in ﬂ uences the transportation of waxy crude oil. The deposited layer grows continuously and hardens during the oil transportation, reducing the e ﬀ ective inside diameter of the crude oil pipeline and the ﬂ ow rate. In extreme cases, the deposited layer may block the crude oil pipeline leading to a loss of production and capital investment. In this paper, wax deposition from multiphase ﬂ ow in ﬁ eld-scale oil pipeline transport systems has been studied. The novelty of this work is to develop a mathematical model that incorporates water-in-oil emulsions, wax precipitation kinetics, molecular di ﬀ usion, and shear dispersion to enable accurate predictions of both the wax deposit growth rate and aging of the deposit. The coupled nonlinear partial di ﬀ erential equations governing the ﬂ ow are discretized in time by a second-order semi-implicit time discretization scheme based on the Adams-Bashforth and Crank-Nicolson methods, which completely decouples the computation of the governing equations. The resulting temporal schemes are discretized in space by the bivariate spectral collocation method based on Chebyshev-Gauss-Lobatto grid points and simulated in MATLAB software to obtain the pro ﬁ les of the ﬂ ow variables. The simulation results are presented in graphical and in tabular forms and discussed. This study found that the deposit thickness is directly proportional to the Reynolds number and inversely proportional to the mass Grashof number, Schmidt number, and Weber number. Deposit aging is rampant during the early stages of wax deposition, after which it stabilizes at a speci ﬁ c value as time elapses. A deposition model to predict the wax deposit thickness and aging is proposed in this study. The ﬁ ndings of this study will help in making informed decisions on the planning of pigging operations, thermal insulation, and other remediation techniques to be applied in controlling wax deposition in ﬁ eld-scale crude oil pipeline systems.


Introduction
A pipeline is the most efficient means for transporting fluids, such as waxy crude oil. The increasing energy demand has caused oil exploration in Africa extends from onshore to offshore fields and to greater depths, such as the Lokichar oil fields in Kenya, where crude oil gets transported through long pipelines in the range of 20-60 km from the reservoir to inland facilities. The long distance for oil transportation means that the time is adequate for oil to cool to temperatures below the wax appearance temperature, usually between 10°C and 30°C. The deposition of solid wax crystals on the inner wall of crude oil pipelines is a challenge in the oil industry. Figure 1 shows a crosssection of a pipeline which is plugged by a significant amount of wax deposit. It presents the general sketch of the problem under study.
Low temperature is the main factor affecting the wax precipitation and deposition process. It means that pipelines in relatively cold places are especially vulnerable to wax deposition. Addressing the problem of wax deposition at an early stage of pipeline design and construction may reduce the operating and maintenance costs of the pipeline system.
Various wax mitigation methods are employed to control wax formation and deposition, such as mechanical techniques (e.g., pigging), heating of the pipeline, and chemical inhibitors. However, heating the crude oil pipeline is not economical in the case of offshore oil exploration. Moreover, excessive usage of additive chemicals can damage the ecosystem in case of accidental release of chemicals to the environment, whereas pigging operation cannot be applied efficiently without prior knowledge of the deposit thickness. In particular, pigs usually get stuck inside the crude oil pipeline due to thick and hard deposits, making the situation worse. In the worst case, production must stop to replace the plugged portion of the crude oil pipeline. This situation is estimated to cost about 40,000,000 dollars per incident [1]. Due to the limitations of the available wax mitigation methods and the persistent challenge of wax deposition, many oil companies are still searching for a better solution to wax deposition.
There are various theoretical and experimental studies in understanding and modeling the formation and deposition of solid wax in crude oil pipeline systems. For example, Singh et al. [2] studied the formation and aging of the gel layer by performing a series of laboratory flow loop experiments. They developed a mathematical model to predict the deposit growth rate and the increase in wax content of the gel with time considering externally cooled pipeline walls. Kok and Saracoglu [3] developed a mathematical model for estimating wax deposition in crude oil pipeline systems. Fusi [4] studied the unsteady flow of waxy crude oil in a laboratory test loop. The flow was analyzed in a nonisothermal condition considering the radial coordinate of the pipe as the only space variable. Nazar et al. [5] experimentally studied the variation of the amount of wax deposition for a mixture of toluene and waxy oil. The study considered both the molecular diffusion of wax as a mechanism of deposition and the sloughing effect due to the hydrodynamic forces of fluid on the deposited wax. Banki et al. [6] studied the mathematical formulation and numerical modeling of wax deposition in crude oil pipelines for laminar flow, using the enthalpy-porosity approach.
Stubsjøen [7] studied analytical and numerical modeling of paraffin wax in pipelines. The study examined the influence of the deposit on the thermal conditions in the pipeline and found that the temperature at the oil-deposit interface increases with an increase in wax deposit thickness. Skjaeraasen et al. [8] studied the pipeline restart behavior at subsea conditions for waxy crude oil, employing advanced computational pressure wave modeling. Fatkhullina et al. [9] studied mathematical modeling of a water-in-oil emulsion droplet behavior under the microwave impact considering unsteady flow. The volume of fluid (VOF) method was adopted to resolve the unsteady dynamics of the free-interface separating the phases, throughout the flow domain.
Zhang et al. [10] developed a wax deposition model to predict the temperature field and location of wax deposition based on a heat-fluid coupling process for a typical rod-tubing-casing-cement-formation with the borehole axis as the center. Ying et al. [11] studied heat transfer analysis of oil phase-change during overhead pipeline shutdown. Magnini and Matar [12] studied wax deposition in crude oil-deposit two-phase flow in a pipeline through interface-resolved numerical simulations. The study assumed that molecular diffusion is the only mechanism responsible for wax deposition and that the effect of water-in-oil emulsions is negligible. The study adopted the volume of fluid (VOF) method to capture the unsteady dynamics of the free oil-deposit interface.
Mouketou and Kolesnikov [13] studied modeling and simulation of multiphase flow applicable to processes in the oil and gas industry. The numerical simulations were performed under steady-state conditions using CFD software. Sun et al. [14] performed an experimental and theoretical study of wax deposition behaviors on a heat-insulated waxy crude oil pipeline in Northeast China. The effects of flow rate and ambient temperature on the thickness and wax content of the deposition layer were investigated. Jiang et al. [15] conducted a numerical study for removing wax deposition by thermal washing for the waxy crude oil pipeline. The study adopted the enthalpy-porosity and volume of fluid (VOF) methods to simulate the melting process of wax in the crude oil pipeline.
From the previous modeling studies on the flow of waxy crude oil in pipeline systems, most of the focus has been put on the single-phase flow of crude oil in pipelines but significant problems may also occur in multiphase flows such as water-oil or gas-oil two-phase flows and water-gas-oil three-phase flows, which are prominent in oil field operations. Water-oil two-phase flows continue to gain interest from the upstream oil industry since the water content of the ore increases significantly with the increase in extraction time. Thus, existing studies have not fully addressed the needs of oil producers since these models have not incorporated the coexistent water during crude oil extraction. The identified research gaps prevent the direct application of existing wax deposition models to field-scale crude oil pipeline systems.
The novelty of this study is to advance the understanding of wax deposition in the presence of water-in-oil emulsions and develop a mathematical model for a wax deposition that incorporates the effects of coexistent water, surface tension force, the porosity of the gel layer, internal heat generation/ absorption, viscous dissipation, precipitation kinetics of  Journal of Applied Mathematics wax, molecular diffusion, and shear dispersion on the growth and aging of the gel layer. The end product is an accurate prediction of both the wax-deposit growth rate and aging (or hardening) of the deposited layer, which is essential for the design of pigging operations to achieve optimum pigging frequency and minimize operational costs. Some specific applications to the present study are the studies by Kadir et al. [16], Anwar Beg et al. [17], and Kuharat et al. [18]. The rest of the paper is organized as follows: Section 2 presents the mathematical model development, Section 3 presents the results of this study, Section 4 presents the validation of the results, Section 5 presents the prediction of wax deposit thickness and aging, and Section 6 presents the summary and conclusions drawn from this study.

Model Formulation
The unsteady two-dimensional flow of waxy crude oil with some water droplets in a model pipeline of circular crosssection, with a semi-infinite length and inner radius R, is considered, as shown in Figure 2. The pipeline is elevated at an angle of α to the horizontal. Cylindrical coordinate system ðr, θ, zÞ is chosen such that r denotes the radial distance from the pipe centerline, θ denotes the tangential direction, and z denotes the stream-wise coordinate or the axial direction. At the time t = 0, waxy crude oil of uniform temperature T ∞ is injected at the pipeline inlet. The pipeline's inner wall is considered smooth, rigid, and impermeable and is maintained at a uniform temperature T wall . The temperature T wall may either be less than or greater than the free-stream (or ambient) temperature T ∞ .
The fluid (i.e., water-oil) is considered a continuous primary phase, whereas solid wax crystals are the discrete secondary phase. Furthermore, the dispersed water droplets flow at the same velocity as the oil phase. The fluid phase and the solid phase are separated with a sharp interface whose outward-drawn unit normal vector is n = ∇ ! ϕ oil /j∇ ! ϕ oil j. The fluid system is described in terms of three pseudocomponents, i.e., oil, wax, and emulsions.
This study employs an interface-tracking technique called the volume of fluid (VOF) model to capture the dynamics of the interface separating the phases. The fluid model adopted makes use of three variables (i.e., ϕ oil , ϕ water , and ϕ gel ) to define the volume fraction of each phase within the flow domain such that The volume fraction of each phase takes values in the interval ½0, 1. If 0 < ϕ i < 1 (for i ∈ foil, water, gelg), then the computational cell contains the interface between the i-th fluid and one or another fluid. If ϕ i = 0, the cell is empty of the i-th fluid. If ϕ i = 1, the cell is full of the i-th fluid.
This study also employs the Pseudosingle Phase (PSP) approach, where the crude oil-water-gel three-phase mixture is treated as a single mixture fluid. The physical properties of the mixture fluid are computed by averaging the corresponding physical properties of the oil, water, and gel and weighted by the corresponding volume fractions, presented in Zheng [19], Yang et al. [20], and Al-Ahmad et al. [21] as where the subscript f denotes "fluid phase" and mix denotes "mixture fluid." Crude oil is a mixture of several hydrocarbon components. Hence, tracking the concentration of every single component is not practical as the flow evolves within the pipeline. Therefore, this study adopts a pragmatic fluid model in which wax exists in two different forms, i.e., dissolved wax and crystallized wax (or precipitated wax). Since the mass concentration of wax (or the amount of wax) which remains soluble (C d ) in the liquid phase must be equal to the total wax concentration (C) minus the amount of wax that has precipitated (C p ), the total mass concentration of wax in the crude oil at any location (r, z) and time t is given by Suppose that the precipitated wax (C p ) exists in two different forms, i.e., aggregated or deposited wax (C a ) and nonaggregated wax (C n ). Thus, the mass concentration of the nonaggregated wax crystals is given by  At temperatures above the WAT, we have Cðr, z, tÞ = C d ðr, z, tÞ since we are dealing with an unsaturated solution. This study focuses on complete saturation; i.e., we suppose that the dissolved wax concentration (C d ) is the one of saturation at all locations in the domain, which depends only on temperature, i.e., To describe the crystalline (or precipitate) component of wax in crude oil, we consider a dimensionless parameter α m known as the aggregation degree of precipitated wax. It represents the ratio of the concentration of aggregated wax to the concentration of precipitated wax [4].
Thus, α m takes values in the interval ½0, 1. It plays a crucial role in determining the specific rheology of the fluid. Eliminating C a between Equations (4) and (6) yields Eliminating C p between Equations (3) and (7), we obtain the relation Equation (8) is used to calculate the concentration of nonaggregated wax crystals (C n ) in the crude oil pipeline once C, C d , and α m are known from the numerical simulations.
The simultaneous action of two wax deposition mechanisms, i.e., molecular diffusion and shear dispersion, is con-sidered. The mass flux of the dissolved wax due to molecular diffusion is given by On the other hand, the mass flux of precipitated wax due to shear dispersion is given by where L 1 and L 2 are nonnegative functions of the shear rate γ n (evaluated for the surface element with normal vector n ! − taken as a unit vector in the radial direction).
The general equations governing the flow of waxy crude oil in pipeline systems are the equations of continuity, momentum, energy, wax concentration, wax precipitation kinetics, oil volume fraction, deposit growth, and deposit aging.
where V ! is the velocity field of the fluid phase, T is the temperature field, R is the inner radius of clean pipe, R eff is the effective radius for oil flow, δ is the thickness of the gel layer, x is the weight fraction of solid wax in the gel layer, and f ðxÞ is the coefficient describing diffusion of wax molecules in the gel layer.
This study assumes that the flux J ! p of the precipitated wax crystals is proportional to the concentration C n of the nonaggregated wax and to the radial velocity gradient ∂u z /∂r only. It is assumed further that the flow is axisymmetric; there is no gas in the crude oil pipeline system, there is no slip of fluid particles at the solid-liquid boundary, and that thermophysical properties in each phase are constants except the variation of the density of the mixture fluid with temperature and species concentration in the body force term so that the usual Boussinesq's approximation is applicable in the boundary layer flow.
Under these assumptions, the specific equations (in cylindrical coordinates) governing the flow of waxy crude oil in field-scale pipeline systems are the equations of continuity, momentum, energy, wax concentration, wax precipitation kinetics, oil volume fraction, deposit growth, and deposit aging.
The corresponding boundary and initial conditions are  To express the governing equations and the initialboundary conditions in their dimensionless forms, the following dimensionless variables are introduced.
Under the transformations (14)- (17), the dimensionless equations governing the flow of waxy crude oil in field-scale pipeline systems are The dimensionless parameters in the model Equations (18)-(26) are defined as follows: α avg = 1:684 − 0:323 ln Q, Journal of Applied Mathematics This study adopts the thermodynamic model for waxy crude oil given in Al-Ahmad et al. [21] as where S f = ð0:0077MW oil − 1:737Þ is the shift factor and M W oil = 6084/ o API − 5:9 is the molecular weight of the crude oil [22]. We consider heavy crude oil whose API gravity is 18°API.
The corresponding dimensionless boundary and initial conditions are The model Equations (18)-(26) together with the initial and boundary Condition (29) are solved numerically using spectral methods. In particular, the discretization is performed in two ways: temporal discretization using a second-order semi-implicit scheme (i.e., a combination of both Adams-Bashforth second-order and Crank-Nicolson) and spatial discretization using bivariate spectral collocation method based on Chebyshev-Gauss-Lobatto points. The resulting numerical schemes are simulated in MATLAB software to obtain the profiles of the flow variables.

Results and Discussion
The results obtained from the present study are presented in graphical and tabular forms and discussed.  Figure 3 shows the radial velocity profiles at different positions. As can be seen, the radial velocity increases as the distance increase from both the pipe inlet and the centerline. The observed trend is because as the flow progresses, the velocity boundary layer thickness reduces, thus accelerating the flow radially. Figure 4 shows the axial velocity profiles at different positions. As can be seen, the axial velocity decreases as the distance increases between the pipe inlet and the centerline. The observed trend is because as the flow progresses, more wax gets deposited on the pipeline's inner wall, which grows and hence reduces the effective flow radius of the pipe. Figure 5 shows the temperature profiles at different positions. As can be seen, the fluid temperature increases as the distance increases from both the pipe inlet and the centerline. The observed trend is because as the flow progresses, the wax deposit absorbs the thermal potential of waxy crude oil due to the high-temperature difference. The outlet temperature begins to decrease and reaches its minimum value. With the reduction of solid wax, the potential of the wax deposit to absorb heat decreases, which increases the outlet temperature. The outlet temperature approaches the inlet temperature and reaches waxy crude oil temperature. Figure 6 shows the total concentration profiles at different positions. As can be seen, the wax concentration increases as the distance increases between the pipe inlet and the centerline. The observed trend is because as the flow progresses, the high-temperature difference leads to thinning of the concentration boundary layer thickness. Hence, more wax molecules diffuse from the bulk of the fluid toward the pipeline wall. Figure 7 shows the wax aggregation degree profiles at different positions. As can be seen, the wax aggregation degree decreases as the distance increases between the pipe inlet and the centerline. The observed trend is because as the flow progresses, more wax molecules precipitate to form solid wax crystals, which aggregate and deposit on the pipeline wall. The increase in the concentration of precipitated wax reduces the aggregation degree. Figure 8 shows the oil volume fraction profiles at different positions. As can be seen, the oil volume fraction decreases as the distance increases between the pipe inlet and the centerline. The observed trend is because the heat gets transferred from waxy crude oil to the wax deposit. The thermal potential of waxy crude oil decreases as it flows.

Journal of Applied Mathematics
Thus, the temperature difference between the wax deposit and waxy crude oil decreases along the pipeline. Therefore, the wax deposit away from the pipe inlet melts slowly at a narrow temperature difference. As the process advances, the temperature of the wax deposit increases, and the tem-perature difference diminishes, which reduces the heat transfer rate. Therefore, the growth rate of the wax deposit reduces.

Deposit Thickness and Aging.
The data (about deposit thickness and weight fraction of wax) extracted from computer simulations are plotted against time while varying the flow parameters. Moreover, the values of the deposit growth and aging rates are presented in Table 1.

Effects of Flow Parameters on Deposit
Thickness. It is observed from Figure 9 that deposit thickness increases with an increase in Reynolds number. The observed trend is due to the enhanced heat convection that promotes large wall temperature gradients. However, large values of the Reynolds number increase the rate of deposit removal. Therefore, the terminal steady-state deposit thickness increases.
It is observed from Figure 10 that deposit thickness decreases with an increase in mass Grashof number. Increasing the mass Grashof number increases the species buoyancy forces in the fluid, which leads to thinning of the hydrodynamic boundary layer and thus accelerating the flow of waxy crude oil in the axial direction. An increase in mass Grashof number increases the deposit removal rate but reduces the    Journal of Applied Mathematics wax deposition rate. Hence, the terminal steady-state deposit thickness decreases as the flow progresses. It is observed from Figure 11 that deposit thickness decreases with an increase in Schmidt number. Increasing the Schmidt number increases the momentum diffusivity, which leads to thinning of the velocity boundary layer and thus accelerating the flow in the axial direction. An increase in Schmidt number increases the deposit removal rate but reduces the wax deposition rate. Hence, the deposit thickness decreases as the flow progresses.
It is observed from Figure 12 that an increase in Weber number results in a decrease in the thickness of the wax deposit. The observed trend is because, at lower surface tension, droplets of deposited phase are continuously detached from the top of the interfacial waves and entrained in the bulk of the flow. At a low Weber number, the droplet entrainment disappears. Hence, the capillary forces promote the formation of a thicker deposit film.

Effects of Flow Parameters on Deposit
Aging. It is observed from Figure 13 that an increase in Reynolds number results in a decrease in the weight fraction of wax in the deposited layer. The observed trend is because the exposure    9 Journal of Applied Mathematics of the deposited layer to a range of temperatures and shear stress triggers physical and chemical processes on the wax deposit, which alters its molecular structure and properties. The thermal gradient across the deposit may result in an internal mass flux which yields a continuous increase in the wax content of the deposited layer. Thus, it leads to the hardening of the deposit as time elapses, a process called aging of the deposit.
It is observed from Figure 14 that an increase in Schmidt number results in a decrease in the weight fraction of wax in the deposited layer. The observed trend is because the Schmidt number decreases the deposit thickness. Hence, the thinner deposit layer increases the temperature gradient in the deposited layer, which increases the diffusion flux in the deposited layer. Thus, it leads to the continuous diffusion of wax molecules into the deposited layer, increasing the wax content in the deposited layer. Table 1 that both the Reynolds number and Schmidt number significantly affect the wax deposit growth and aging rates. As the wax flow progresses, the deposit growth and aging rates decrease due to an increase in the interfacial temperature (because the oil-deposit interface moves away from the pipeline wall) and to the onset of deposit removal.

Validation of the Results
The results of this study are validated by considering the wax deposition data of Magnini and Matar [12]. The authors studied wax deposition in crude oil flows in a pipeline via interface-resolved numerical simulations. A comparison of the results on deposit thickness profiles in Magnini and Matar [12] with those obtained from the present study is displayed in Figures 15(a) and 15(b), respectively.
The results reveal the same trend in the deposit thickness as time elapses; i.e., from the onset of the flow, the total amount of deposit within the pipeline increases until deposition and removal rates balance out and the deposit thickness reaches a steady state. The observed trend is due to the governing effect of heat convection at the early stages of wax deposition, where the radial temperature gradient increases with the flow rate. At later stages of wax deposition, i.e., when the deposit is thicker, higher flow rates promote shear-induced deposit removal, thus yielding lower deposit thicknesses.

Prediction of Wax Deposition Thickness and Aging
This study proposes simplified prediction models for the deposit thickness and aging by the Lagrange interpolation technique. Using MATLAB toolkit for basic fitting applied to Figures 15(b) and 13, we obtain the following Lagrange polynomial of degree 4 for estimating the deposit thickness and aging as functions of time t.

Summary and Conclusions
A mathematical model has been developed to predict wax deposition and aging from multiphase flow in field-scale crude oil pipeline transport systems using both the VOF and PSP approaches. Waxy crude oil and the wax deposit are two immiscible phases separated by an interface. The wax concentration has been modeled by considering a pragmatic fluid model in which wax exists in two different forms, i.e., dissolved wax and crystallized wax (or precipitated wax). Molecular diffusion and shear dispersion have been considered the only deposition mechanisms. The model is a system of coupled nonlinear partial differential equations, i.e., the continuity, momentum, energy, species transport, wax aggregation degree, conservation of oil volume fraction, deposit growth, and deposit aging equations, for the wax-fluid two-phase flow in the presence of waterin-oil emulsions. The model equations are solved numeri-cally by the bivariate spectral collocation method. From the study, the following conclusions have been drawn: (i) Initially (i.e., at time t = 0), no deposit is present, and thus removal does not occur yet. At the onset of the flow, wall deposition starts due to the radial temperature gradient and because the wall temperature is below the WAT. As time elapses, the deposited layer grows because the wax deposition rate exceeds the deposit removal rate. When the deposit is sufficiently thick such that deposition and removal rates balance out, the average thickness of the deposit becomes constant. The steady-state values of the deposition rates depend strongly on the film thickness, i.e., the balance between deposition and removal rates, which determines the temperature of the crude-deposit interface  Figure 15: Graphs of deposit thickness in Magnini and Matar [12] and in the present study.