Reactive Transport and Its Implications on Heavy Oil HTGC Analysis – A Coupled Thermo-Hydro-Chemical (THC) Multiphysics Modelling Approach

Abstract

This chapter provides an insight into the reactive transport in a capillary column which heavy-oil hydrocarbons undergo when analysed by high temperature gas chromatography (HTGC), and their implications on characterisation outcomes, namely thermal cracking of the injected sample; and incomplete or non-elution of heavy components from the column, by using a coupled Thermo-Hydro-Chemical (THC) multiphysics modelling approach. For this purpose, a computational coupled THC, multicomponent, multi-physics model is developed, accounting for: multiphase equilibrium using an in-house, extended thermodynamics distribution factors dataset, up to nC98H198; transport and fluid flow in COMSOL and MATLAB; and chemical reactions using kinetics and mechanisms of the thermal cracking, in CHEMKIN. The determination of the former extended dataset is presented using two complementary HTGC modes: i) High-Efficiency mode, with a long column operated at low flow rate; and ii) true SimDist mode, with a short column operated at high flow rate and elution up to nC100H202.

1. Introduction

The objective of this chapter is to understand the reactive transport occurring during the High Temperature Gas Chromatography (HTGC) analysis of heavy oil hydrocarbons at common conditions and thereby quantify the implications on the final characterisation results in terms of: (i) the degree of thermal cracking of the original sample; and (ii) the non-elution of heavy components from the HTGC column by using a combined Thermo-Hydro-Chemical (THC) coupled multiphysics modelling approach.

For this endeavour, a synergy between experimental and computational coupled multi-physics approaches, has been carried out to account for three physicochemical processes: thermodynamic equilibrium fluid-flow; transport of chemical species; and chemical reactions.

An outline is given of the experimental approach used, with explanation of the methodology for extending the distribution factors data-set, necessary to describe the first process.

On the computational side, an in-house coupled multi-physics model has been developed using MATLAB [1] as language host, to couple the above three processes. The former is described, using as input to the multi-physics model the extended, experimental distribution factors dataset: and the latter two processes are described using: COMSOL Multi-physics [2] and MATLAB, and CHEMKIN [3] respectively.

Finally, the implication of the inter-related, multi-physics, physicochemical processes is discussed, based on the results from the coupled THC multi-physics model.

2. Experimental overview

2.1 Outline of HTGC methodology

Detailed accounts of the experimental procedures have been published previously [4], covering the generation of isothermal and temperature-programmed retention data for n-alkanes and polywax mixtures, on poly-dimethyl-siloxane HTGC columns, up to 430°C. (i.e. based on well-established SimDist techniques). This database then enabled the distribution factors of heavy n-alkanes to be derived up to nC₉₈H₁₉₈, which were unavailable in the literature.

2.2 Methodology for extending distribution factors up to nC₉₈H₁₉₈

In the absence of a database of distribution factors of heavy n-alkanes, it was necessary to obtain insight into their behaviour inside the HTGC column, requiring development of a comprehensive methodology to extend the existing, limited amount of data up to around nC₉₈H₁₉₆ [4].

Hernandez et al. [4] derived a temperature-dependent function of distribution factors which has been applied to a series of n-alkanes spanning (nC₁₂H₂₆ - nC₉₈H₁₉₆) by combining Eq. (1) and numerous isothermal experiments carried out using two SGE HT5 GC capillary columns [5] of different lengths, and under two HTGC methods as follows:

1. For the long column, low flow rates were used, to obtain efficient resolution of eluting n-alkanes, spanning the range (nC₁₂H₂₆–nC₆₄H₁₃₀), under constant inlet pressure measurements conditions.

2. For the short column, ASTM method D7169–11 [6] was applied to achieve extended SimDist analysis, spanning the range nC₁₂H₂₆–nC₉₈H₁₉₈ under constant flow rate measurement conditions.

In both columns, the standard samples (ASTM D5442) was used and at least 3 isothermal GC measurements were carried out from 80 to 420°C at 20°C intervals, and lastly at 430°C. Further details can be found in [4].

3. Theory

3.1 Physicochemical HTGC workflow

In order to understand the Reactive Transport (inter-related Thermo-Hydro-Chemical multi-physics processes) occurring during HTGC analysis of heavy oils, it is necessary to consider them step-wise, within the column.

• as the vaporised hydrocarbon/CS₂ solvent sample is transported through the column, the Diffusion Process of each component is considered to be negligible in comparison with the Convection Process [7].

• the stationary phase controls the adsorption of each component of the sample as a function of its boiling point in relation to the oven temperature program

• the vaporised hydrocarbon sample encounters the stationary phase, where each component experiences a gas–liquid thermodynamic equilibrium process at the prevailing temperature and pressure, between the stationary and mobile phases.

• Each component of the sample which has been dissolved in the stationary phase is retained at the same point until its vapour pressure equals that of the carrier gas, as the column temperature program proceeds, and is then transported to the GC outlet. It is assumed that no diffusivity effects occur in the liquid phase.

• The solvent has a low boiling point compared to the hydrocarbon mixture, being the first component to be eluted or released by the stationary phase and transported by the carrier gas to the column outlet. It therefore generates the first peak in the gas chromatogram results, representing the least retained solute in the sample.

• Each component might encounter other components in the carrier gas and it is assumed that all of them are transported as a batch mixture.

• Each mixture is at risk of suffering a thermal cracking process at a given temperature in the gas phase, depending on whether the kinetics and mechanism at a given temperature reach the necessary energy to break some of the chemical bonds of the component, creating smaller chemical species.

• Each component of the sample whose boiling point is not reached during the temperature program, is retained by the stationary phase, and at risk of either incomplete elution, or non-elution from the column.

• The longer a heavy component is retained by the stationary phase before desorbing, the greater the risk of thermal cracking due to its greater exposure to the highest temperature.

3.2 Thermo-hydro-chemical processes

Thus, the three physicochemical processes in a heavy-oil HTGC analysis considered in this work are:

• Multiphase equilibrium:

  Solvation thermodynamics using experimental data.

• Transport and fluid flow:

  Convection in MATLAB

• Chemical reactions of thermal cracking:

Kinetics and mechanisms of simplified mixture in CHEMKIN,

which fall within the classification Thermo-Hydro-Chemical multi-physics [8].

3.2.1 Multiphase equilibrium: solvation thermodynamics

The basis of the gas chromatography separation process centres on the non-isothermal multiphase equilibrium of each of the chemical species in the mixture sample between the stationary phase and the gas phase (transported by the carrier gas) taking advantage of their distinct boiling points.

This equilibrium is established in multiple stages throughout the length of the capillary column. The analysis mixture sample is dissolved and retained in the stationary phase at low initial temperatures and each component comprising the mixture is evaporated and separated from the sample mixture once its boiling point temperatures and pressure is reached. Thus, solvation thermodynamics is used to describe the gas–liquid equilibrium of each chemical species inside the GC column.

The temperature-dependent expression for the distribution factor, K, is described in Eq. (1) and was obtained by solving the thermodynamic equilibrium of the solvation of a solute in the bulk solvent [9] expressed in terms of the Gibbs free energy at a given temperature and by the logarithm of the solute molecule’s numeral density ratio in both phases [10, 11] or the ratio between the molar concentration of the two phases.

The mass transfer is assumed to be governed only by the interaction between the solute and the stationary phase, while the interactions between solute-solute and solute-carrier gas are assumed to be negligible as the interfacial and extra-column effects leading to non-equilibrium conditions [12].

A semi-empirical model [13, 14] developed by Castells et al. [15] for the determination of the isothermal retention times as function of the hold-up time, t_M and the solvation time expressed by the Gibbs free energy is expressed in the terms of ΔH and ΔS, which represent the changes in enthalpy and entropy associated with the transfer of solute from the stationary phase to the mobile phase at a given temperature T.

$$K\left(T\left(t\right)\right)=\beta \left[\frac{t_r}{t_M}-1\right]=exp\left[a_0+a_1\frac{1}{T\left(t\right)}\right]$$

$$a_0=\frac{\mathit{\Delta S}\left(T\right)}{R};a_1=-\left(\frac{\mathit{\Delta H}\left(T\right)}{R}\right)$$

$$\beta =\frac{{\left(2r_o-2w\right)}^2}{2{r_o}^2-{\left(2r_o-2w\right)}^2}$$E1

In Eq. (1), K corresponds to the distribution factor, with t_r and t_m representing the retention time of the solute and the hold-up time, respectively. β is the phase ratio of the column, r_o and w correspond to the inner radius of the column, and the film thickness of the stationary phase. ΔH and ΔS, correspond to the delta changes in enthalpy and entropy associated with the transfer of solute from the stationary phase to the mobile phase.

Aldaeus [16] has proposed two retention mechanisms according to the nature of the separation hold between the analyte and the stationary phase, based on the semi-empirical values of the thermodynamic properties of Eq. (1).

3.2.2 Transport and fluid flow: convection

The Snijders [17] method for calculating the retention times in gas chromatography is based on the peak position determination which is not affected by the diffusion effects but by the convection process only [16].

Eq. (2) expresses the convection process in terms of the effective velocity v_eff of the analyte in the carrier gas. Discretized into finite time-steps of Eq. (2) allows tracking of the position of the analyte at every x position through the GC column at every time step until the peak reaches the column outlet [18, 19] at the final time step which cumulated represents the retention time for that analyte as explained in [4]. And lastly, K and β are the distribution factor and phase ratio of the column described in Eq. (1) and v_m is the velocity of the mobile phase.

$${\upsilon }_{\mathit{eff},\mathit{i}}\left(x,t\right)=\frac{{\upsilon }_M\left(x,t\right)}{1+\frac{K_i\left(T\right)}{\beta }}$$E2

Integrating the differential form of the Hagen-Poiseuille fluid mechanics Equations [10, 18] through the length of the column allows calculation of v_m as described in Eq. (3). This expression relates the carrier gas velocity to the pressure gradient at any position in the column [18] by a proportional constant. The latter depends of the geometry of the cross-section which in this case is for a column of circular cross-section [20]:

$${\upsilon }_M\left(x,t\right)=\frac{{r_o}^2·\left({P_{\mathit{\text{in}}}}^2-{P_{\mathit{\text{out}}}}^2\right)}{16·\eta \left(T\left(t\right)\right)·L·P\left(x\right)}$$E3

In Eq. (3), η(T(t)) corresponds to the viscosity of the carrier gas [21, 22]. (See summarised details in [19]), P_in and P_out are the inlet and outlet pressures of the GC column. r_o is the inner radius of the column and P(x) is the pressure at position x described with Eq. (4), which is obtained by integrating the Hagen-Poiseuille equation between the inlet and outlet position, of a differential element and assuming incompressibility of the gas in each element at position x, due to the extremely low pressure-drop in gas chromatography [10].

$$P\left(x\right)=\sqrt[2]{{P_{\mathit{\text{in}}}}^2-\left({P_{\mathit{\text{in}}}}^2-{P_{\mathit{\text{out}}}}^2\right)\frac{x}{L}}$$E4

3.2.3 Chemical reactions: kinetics and mechanism of thermal cracking

The large number of species in the reduced free-radical pyrolysis model developed in [23] has imposed a need to develop a reduced molecular pyrolysis model, comprising 11 n-alkanes (nC₁₄H₃₀, nC₁₆H₃₂, nC₂₀H₄₂, nC₂₅H₅₂, nC₃₀H₆₂, nC₃₅H₇₂, nC₄₀H₈₂, nC₅₀H₁₀₂, nC₆₀H₁₂₂, nC₇₀H₁₄₂, and nC₈₀H₁₆₂).

In this work, a “class” molecular mechanism has been obtained after applying the following three rearrangements to our reduced molecular mechanism model [7]:

• Lumping of molecules belonging to the global class “C₁₅ plus” which are produced by an n-alkane reactant.

• Lumping of n-alkane reactants, which produced n-alkane reactants or lighter class.

• Lumping of global class C₁₅ as reactant.

Refer to [23], to understand the thermal cracking original kinetic and mechanism model development, and refer to [7] to understand the detailed explanation of the kinetics and mechanism reduction procedure from molecular mechanism model to a “class” molecular mechanism.

The optimised reduced “class” molecular mechanism used in this work is composed of 127 molecular reactions and 17 species (11 n-alkanes, and 6 “class” molecular pyrolysis products) which has been obtained after applying to the whole mechanism the above rearrangement and its corresponding kinetic data [7].

Thus, the final reduced molecular mechanism, accounts for:

• 11 original n-alkanes (reactants): nC₁₄, nC₁₆, nC₂₀, nC₂₅, nC₃₀, nC₃₅, nC₄₀, nC₅₀, nC₆₀, nC₇₀, nC_80.

• 6 classes: alkene, CH₄, C₂H₆, C₃-C₅, C₆-C₁₃ and C₁₅ plus.

Finally, as summarised in Table 1 the number of reactions of the original free-radical pyrolysis model has been reduced from 7055 to 127, and the number of species from 336 to 17, whilst still yielding very good accuracy as depicted in Figure 1.

	Pyrolysis mechanisms
	Radicals	“Class” Molecular
Mixture of Heavy Oils	nC14, nC16, nC20, nC25, nC30, nC35, nC40, nC45, nC50, nC55, nC60, nC65, nC70, nC75, nC80	nC14, nC16, nC20, nC25, nC30, nC35, nC40, nC50, nC60, nC70, nC80
Reactants	15	11
Reactions	7055	127
Species	336	17
Molecules	242
Radicals	94

Table 1.

Summary of size of the mechanistic kinetics models developed.

Figure 1.

Comparison of free radical model and “class” molecular model for heavy n-alkanes mixtures. (simulation of a closed reactor at 1 MPa).

4. Computational multi-physics

An understanding of the Thermo-Hydro-Chemical (THC) processes occurring inside an HTGC column during the analysis of heavy oil hydrocarbons was obtained through detailed study with an in-house coupled THC model.

The coupling of the physico-chemical processes is sequential due to the complexity of the system, and the level of detail with which each process has been described. Hence, a fully coupled model is prohibited while a sequential coupling can handle effectively the following processes:

• the thermodynamics equilibrium is described using experimental input data of a series of n-alkanes spanning (nC₁₂H₂₆-nC₉₈H₁₉₆) [4].

• the chemical reactions are described using kinetics and mechanistic modelling of 11 n-alkanes: nC₁₄, nC₁₆, nC₂₀, nC₂₅, nC₃₀, nC₃₅, nC₄₀, nC₅₀, nC₆₀, nC₇₀, nC₈₀ and 6 class molecules: alkene, CH₄, C₂H₆, C₃-C₅, C₆-C₁₃ and C₁₅ plus.

• the convection process is described using the Hagen-Poiseuille fluid mechanics equations [10, 18].

The sequence of these processes was arranged using short time intervals where the temperature was constant during the temperature ramp, and with a batch size as described using Golay’s theory [24] for diffusion and convection processes.

4.1 Computational HTGC workflow

From a computational modelling perspective, the multi-physics processes can be simplified and described as follows:

1. The system to model is delimited to the gas phase inside the HTGC column, comprising only the carrier gas transporting each component from the column inlet to the outlet.

2. The control volume is the inner volume of the coiled capillary GC column (e.g.,typically 0.53 mm internal diameter, and 25-30 m length).

3. The fluid flow and transport of chemical species considers only the convection process of the carrier gas transporting each of the species at its own speed from the GC inlet to the GC outlet. The diffusion of the chemical species in the gas phase has been concluded to be negligible based on a previous investigation [19] which compares the advection process (convection plus diffusion) and the convection process only, demonstrating no difference between the chromatogram peaks results using both approaches.

4. The equilibrium occurring in the interface between the stationary phase and the gas phase is modelled through the extended thermodynamics distribution factors dataset obtained previously [4], for each one of the heavy n-alkanes hydrocarbons mixture.

5. The heavy oil mixture is simplified to one comprising only heavy n-alkanes, ranging from nC₂₅H₅₂ to nC₉₈H₁₉₈ which is suitably representative, bearing in mind that long-chain n-alkanes are the most susceptible to thermal cracking.

6. Their thermal cracking is modelled in the gas phase only, using a reduced and optimised kinetics and mechanism developed previously [7] and validated against a detailed mechanism of the heavy oil mixture developed initially [23].

4.2 Coupled multi-physics workflow

Finally, the coupling of the three physics involved is made in a sequential order as follows:

1. The column is treated as a series of batches of 1 mm

2. The time for the multi-physics processes to occur in each batch is the time for each species to travel at the carrier gas velocity, from the centroid of one batch to the following one, with velocity derived from volumetric flow measurement and column i.d.

3. The chemical species are modelled through convection only, using MATLAB. No diffusivity is included in this final work as explained in [7]

4. The chemical species undergo an equilibrium which is reached based on the thermodynamics of solvation using the extended distribution factors data set introduced previously [4]

5. The position of each of the chemical species is calculated and each batch is created according to the chemical species found in the batch volume. Thus, the chemical reactions of the thermal cracking are modelled in each mini-batch reactor.

6. Each batch is updated with the new species found inside.

7. The convection of the new chemical species is modelled as described in c) and so on, until all the batches arrive at the GC outlet sequentially.

4.3 Discretization methods

This work uses the discretization method introduced by Snijders [17], which predicts the peak width of the solute zone as the space that a solute migrating through the column occupies [25]. This approach of the convection model was successfully coupled to the reduced molecular pyrolysis model from [7].

Equal time segments are used to discretize the simulation as proposed by Snijders [17] for enabling isothermal properties to be used at every time-step according to the ramp of temperature used. Thus, a sufficiently small time-step permits a uniform pressure to be assumed in the column segments traversed by the solute.

The local plate height (H) is calculated at every time-step based on the Golay [24] expression for open tubular columns, as shown in Eq. (5), where k is the retention factor.

$$\begin{array}{c}H\left(x,t\right)=2·\frac{D_M\left(x,t\right)}{{\upsilon }_M\left(x,t\right)}+{\upsilon }_M\left(x,t\right)\left\{\left[\frac{1+6·k\left(T\left(t\right)\right)+11·k{\left(T\left(t\right)\right)}^2}{24·{\left[1+k\left(T\left(t\right)\right)\right]}^2}·\frac{{r_o}^2}{D_M\left(x,t\right)}\right]\right.\\ \left.+\left[\frac{2·k\left(T\left(t\right)\right)}{3·{\left[1+k\left(T\left(t\right)\right)\right]}^2}\frac{w^2}{D_s\left(x,t\right)}\right]\right\}\end{array}$$E5

Note here that k is dimensionless, being derived from the distribution factor, K, and the phase ratio of the column, β namely K/β, with K corresponding to the ratio between the (moles/volume) in stationary phase to the (moles/volume) in gas phase. r_o and w correspond to the inner radius of the column and the film thickness of the stationary phase D_s, and D_m correspond to the diffusion constant respectively in the stationary and mobile phase. v_m corresponds to the velocity of migration of the carrier gas.

The local zone variance (σx²) in the distance of a solute from the zone centroid at a given position x, can be calculated using Eq. (6), representing the solute band’s spreading.

$${{\sigma }_x}^2\left(\Delta x_n\right)=H\left(x_n,t_n\right)·\Delta x_n$$E6

Eq. (7) describes the increment in the zone variance length, and can be obtained by the summation of all the local contributions of zone variances. Giddings [26] explained that at every time step, the correction is applied for the expansion of the solute zone due to the reduction in pressure (P) along the column.

$${{\sigma }_x}^2\left(x_n\right)=\left[\sum _{i=1}^{n-1}{{\sigma }_x}^2\left(\Delta x_i\right)\right]·\left\{\frac{P\left(x_{n-1}\right)}{P\left(x_n\right)}\right\}+{{\sigma }_x}^2\left(\Delta x_n\right)$$E7

This approach, has been programmed in MATLAB, and has been compared in [7] with the solution yielded by the COMSOL-MATLAB model developed in [19], which solves the diffusive-convective equation by finite elements.

The comparison study confirmed excellent agreement in predictions of the zone’s centroid (average relative error of 1.1%) and of the zone’s standard deviations (average relative error of 3%), as depicted in Figure 2.

Figure 2.

Comparison of zone standard deviation and zone centroid of nC₁₂H₂₆, predicted using an iterative analytic approach¹¹ using MATLAB and solving the diffusive-convection equation by finite element in COMSOL. (Column dimensions Table 3 and temperature programming Table 2).

Thus, the analytical model implemented in MATLAB has been coupled to the reduced molecular pyrolysis model described above, and as detailed in [7], by calling CHEMKIN at every time step iteration, and using feedback between the two models until each component elutes from the GC column.

4.4 Coupled thermo-hydro-chemical processes

Both the reduced molecular pyrolysis model and the analytic iterative GC model derive from the prior high-performance improvement process required for ultimately attaining an efficient coupled physics-chemical model.

The latter can predict the zone’s centroid, the standard deviation and the pyrolysis decomposition of every solute studied for both as a mixture and as a single component according to the position of every solute related to the batch width at every time-step.

In order to maintain a constant temperature at every time-step, a constant time-step has been implemented, permitting an increment of 1°C every 4 seconds (corresponding to the ramp of 15°C/min, used).

Initially, for every component studied, the position of the zone’s centroid in the next time step (x_i+1), is calculated, using Snijders [23, 27] approach Eq. (8) (see ref. [19]), the distribution factor (K), and the phase ratio (β).

$$x_{i+1}=x_i+\frac{{\upsilon }_M\left(x_i,t_i\right)}{1+\frac{K_i\left(T\left(t_i\right)\right)}{\beta }}·\mathit{\Delta t}$$E8

Figure 3, shows the algorithm explaining the global calculation carried out by the coupled THC model for an heavy oil analysis by HTGC, using the above models as explained previously. For more detail refer to [7].

Figure 3.

Algorithm of the pyrolysis-GC coupled model.

5. Implications: results

The implications of the THC processes during an HTGC heavy oil analysis can be summarised under two headings:

Thermal cracking risk of the original sample.

Non-elution or incomplete elution of the sample.

A detailed analysis of these implications is presented, based on the results of the in-house developed THC multiphysics model, described above.

5.1 Thermal cracking risk of the original sample

The cumulative conversion due to pyrolysis of the 11 n-alkanes is studied in [4], in order to analyse their risk, as depicted in Figure 4. For each component the ratio is calculated of the cumulative mass lost (Figure 5) due to thermal cracking, compared to the mass injected.

Figure 4.

Accumulative mass lost due to thermal cracking for n-alkanes (nC₁₄, nC₁₆, nC₂₀, nC₂₅, nC₃₀, nC₃₅, nC₄₀, nC₅₀, nC₆₀, nC₇₀, nC₈₀) at a common HTGC temperature programming (Table 2) in a HT5 column with dimension summarised in (Table 3).

Figure 5.

Cumulative conversion due to thermal cracking for n-alkanes (nC₁₄, nC₁₆, nC₂₀, nC₂₅, nC₃₀, nC₃₅, nC₄₀, nC₅₀, nC₆₀, nC₇₀, nC₈₀) at a common HTGC temperature programming (Table 2) in a HT5 column with dimension summarised in (Table 3).

As would be expected, no pyrolysis reaction occurs in the case of nC₁₄H₃₀ and nC₁₆H₃₄ with the temperature program used (Table 2), and their associated short residence times inside the GC column. Similarly, within the range nC₂₀H₄₂ to nC₄₀H₈₂, insignificant conversion occurs, whereas in the case of nC₅₀H₁₀₂ the maximum mass loss through thermal decomposition before elution is 0.003%.

Parameters	Values
Tinitial (°C)	10
Hold up time at Tinicial (min)	0
ramp of T (°C/min)	15
Tmax (°C)	430
Hold up time at Tmax (min)	12

Table 2.

Temperature programming.

	SGE HT5 GC Column
Length [m]	12
Diameter [mm]	0.53
Film thickness [um]	0.15

Table 3.

Column dimensions of in-house HTGC.

Low but detectable mass loss occurs with the heaviest n-alkanes. nC₆₀H₁₂₂ has a significant loss in the stationary phase where only 2.43·10⁻¹² g are released to the gas phase with the remainder trapped in the stationary phase. Further, pyrolysis loss begins at 373°C with a 0.001% cumulative mass conversion. nC₇₀H_142, presents a cumulative conversion of 0.001% at 385°C with only 2.32·10⁻¹⁰ g released in the gas phase and the rest trapped in the stationary phase.

It should be noted that at the time-step when nC₆₀H₁₂₂ decomposition starts, nC₅₀H₁₀₂ is virtually totally eluted (99.9%) eluted, and hence the pyrolysis products present no risk of co-elution with the latter. Rather, the pyrolysis products of nC₆₀H₁₂₂ are released gradually, evidenced by a slowly increasing baseline signal.

Similarly, nC₇₀H₁₄₂ starts to decompose when located 1.02 m away from the GC inlet, and 0.68 minutes after nC₆₀H₁₂₂ is essentially fully eluted (99.99%). Therefore, the pyrolysis products present no risk of co-elution with, nor distortion of the peak for nC₆₀H₁₂₂.

Lastly, nC₈₀H₁₆₂ starts to decompose at 0.41 m from the column inlet, while nC₇₀H₁₄₂ is located 1.64 m from the inlet. Thus, when nC₇₀H₁₄₂ is essentially fully eluted (99.99%), at 7.83 m from the column inlet, nC₈₀H₁₆₂ has undergone a cumulative conversion of 0.52% mass loss by pyrolysis, relative to mass injected. That equates to 3.97.10⁻¹¹ g of nC₈₀H₁₆₂ converted into pyrolysis products, and which co-elutes with nC₇₀H₁₄₂, resulting in unreliable quantification.

5.2 Non-elution of heavy components from the column

For the determination of non/incomplete elution of heavy n-alkanes, the data set of distribution factors of the n-alkanes spanning the range from nC₁₂H₂₆ to nC₉₈H₁₉₈, [4] was used as main input for the calculation of the degree of elution of each of the n-alkanes studied.

The degree of elution has been introduced in order to determine the non/incomplete elution of heavy n-alkanes (as explained in [19]) as depicted in Figure 6.

Figure 6.

Degree of elution vs. transit time of each component “i”: n-alkanes in the range of C₁₄H₃₀ to nC₈₀H₁₆₂. Degree of elution = moles of “i” inside the GC column at time (t) /moles injected of “i”.

Alkanes heavier than nC₆₀H₁₂₂ elute during the isothermal plateau of the temperature programmed at 430°C. Therefore, constant distribution factors apply for the re-equilibration period, when characteristic peak broadening is observable. (c.f. the essentially symmetrical peaks associated with temperature programmed GC analyses).

nC₇₀H₁₄₂ starts to elute at 29 minutes with a 99.99% degree of elution at 31.3 minutes, and 100% at 31.5 minutes. nC₇₀H₁₄₂ takes 2.5 minutes to elute completely.

nC₈₀H₁₆₂ starts to elute at 33.8 minutes, with a degree of elution of 99.99% at 40.9 minutes and 100% after 42.3 minutes. nC₈₀H₁₆₂ takes 7.1 minutes to elute and 8.5 minutes to completely elute.

In this simulated study, components from nC₇₀H₁₄₂ and above, elute so slowly that peak resolution for the group cannot be assessed. Rather, in practice, a continuum is observed, in the form of a gradually increasing baseline, rising to a plateau which gradually reduces during the final isothermal period of the oven temperature program.

It is interesting to note that 99.99% of nC₈₀H₁₆₂ requires to elute 12.9 minutes at the isothermal conditions at the maximum temperature (430°C) of the analysis. of 99.99%. This means that this component is not normally taken into account in the GC calculations, due to the shorter period of time and stationary phase bleeding.

6. Conclusions

This chapter provides an insight into the analysis of the Reactive Transport process occurring during the analysis of heavy oil hydrocarbons inside a High Temperature Gas Chromatography column, and the implication that those interrelated physicochemical processes generate, by application of a Thermo-Hydro-Chemical (THC) coupled multiphysics approach.

The number of species in the reduced free-radical pyrolysis model developed in [19] has imposed a need to develop a reduced molecular pyrolysis model, comprising 11 n-alkanes (nC₁₄H₃₀, nC₁₆H₃₂, nC₂₀H₄₂, nC₂₅H₅₂, nC₃₀H₆₂, nC₃₅H₇₂, nC₄₀H₈₂, nC₅₀H₁₀₂, nC₆₀H₁₂₂, nC₇₀H₁₄₂, and nC₈₀H₁₆₂). The number of reactions has been reduced from 7055 to 127, and the number of species from 336 to 17, whilst still yielding very good accuracy.

THC multi-physics model has been implemented to resolve the HTGC limitations. The cumulative pyrolysis conversion of the 11 n-alkanes studied in this work, suggests that 0.52% of the mass injected of nC₈₀H₁₆₂, thermally decomposed before nC₇₀H₁₄₂. Therefore, co-elution of nC₇₀H₁₄₂ and the pyrolysis product of nC₈₀H₁₆₂ makes the GC analysis of nC₇₀H₁₄₂ and heavier n-alkanes no longer reliable.

The degree of elution of the 11 n-alkanes studied in the chapter confirms that alkanes heavier than nC₇₀H₁₄₂ take progressively longer to elute completely from the column, viz. nC₇₀H₁₄₂ takes 2.3 minutes and nC₈₀H₁₆₂ takes 7.1 minutes, with co-elution of decomposition products in each case compromising their analyses.

Finally, nC₈₀H₁₆₂ takes 12.9 minutes to completely elute during the isothermal plateau, resulting in no distinct peak being observable. Consequently, the eluting component will be masked in the FID plateau signal, in combination with column bleed products. As a result the nC₈₀H₁₆₂ analysis may not be utilised under these HTGC conditions.

Acknowledgments

The authors wish to thank the members of our JIP: Marathon Oil Corporation, Schlumberger and Total for both their technical and financial support during this project.

Deff	Effective average Diffusivity (unit length2/unit time)
D	Apparent diffusion coefficient, represents all factors causing dispersion (unit length2/unit time)
DM	Diffusion constant, mobile phase (unit length2/unit time)
DS	Diffusion constant, stationary phase (unit length2/unit time)
H(x,t))	Column plate height, spatial rate of dispersion of a zone (unit length)
K	Distribution factor of a compound (moles/volume) in stationary phase/(moles/volume) in gas phase)
k	retention factor of a compound (moles in stationary phase/moles in gas phase).
L	Length of the GC column (unit length)
m(x,to)	mass profile for every analyte (particles/unit length).
Ni,M	Moles of component “i” in the mobile phase.
Ni,S	Moles of component “i” in the stationary phase.
P(x)	Pressure at position x (Pa)
Pin	Pressure at the GC colum inlet.
Pout	Pressure at the GC column outlet (Pa)
ro	Internal radius of GC column (unit length)
ramp T	Ramp of temperature of the temperature programmed.
T(t)	Temperature at the time t.
To	Initial temperature of the temperature programmed.
t	time (unit time)
veff	Effective cross-sectional average velocity (unit length/unit time)
vM	Velocity of migration of the carrier gas (unit length/unit time)
w	Film thickness (unit length)
Xi	Fraction of component i in the gas phase relative to the moles in both stationary and gas phase.
x0	Centroid of Gaussian distribution of distribution of component inside the GC column (unit length)
x	Position of the component’s dispersal around the centroid x0.(unit length)
Greek letters
σ	Standard deviation of the distribution of component inside the GC column (unit length)
β	Phase ratio (volume of mobile phase in the column to the volume of stationary phase).
ηm	Viscosity of the carrier gas.(μPa·s).
Δt	Time step (unit time).