Open access peer-reviewed chapter

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

Written By

Diana Margarita Hernandez-Baez, Alastair Reid, Antonin Chapoy, Bahman Tohidi, Roda Bounaceur and François Montel

Submitted: January 22nd, 2021 Reviewed: May 28th, 2021 Published: August 5th, 2021

DOI: 10.5772/intechopen.98614

Chapter metrics overview

96 Chapter Downloads

View Full Metrics


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.


  • Reactive Transport
  • Thermo-Hydro-Chemical (THC) processes
  • coupled THC modelling
  • coupled multi-physics
  • multiphase equilibrium
  • solvation thermodynamics
  • transport and fluid flow
  • chemical reactions
  • kinetics and mechanism of thermal cracking
  • pyrolysis
  • Heavy n-alkanes thermodynamics distribution factors
  • High-temperature gas chromatography (HTGC)
  • heavy-oil characterisation
  • Gas Chromatography modelling
  • coupled multiphysics modelling

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 nC98H198, which were unavailable in the literature.

2.2 Methodology for extending distribution factors up to nC98H198

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 nC98H196 [4].

Hernandez et al. [4] derived a temperature-dependent function of distribution factors which has been applied to a series of n-alkanes spanning (nC12H26 - nC98H196) by combining Eq. (1) and numerous isothermal experiments carried out using two SGEHT5 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 (nC12H26–nC64H130), 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 nC12H26–nC98H198 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/CS2 solvent sample is transported through the column, the Diffusion Processof 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 processat 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 processat 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, tMand the solvation time expressed by the Gibbs free energy is expressed in the terms of ΔHand Δ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.


In Eq. (1), Kcorresponds to the distribution factor, with trand tmrepresenting the retention time of the solute and the hold-up time, respectively. βis the phase ratio of the column, roand wcorrespond to the inner radius of the column, and the film thickness of the stationary phase. ΔHand Δ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 veff 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 xposition 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, Kand βare the distribution factor and phase ratio of the column described in Eq. (1) and vmis the velocity of the mobile phase.


Integrating the differential form of the Hagen-Poiseuille fluid mechanics Equations [10, 18] through the length of the column allows calculation of vmas 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]:


In Eq. (3), η(T(t))corresponds to the viscosity of the carrier gas [21, 22]. (See summarised details in [19]), Pinand Poutare the inlet and outlet pressures of the GC column. rois the inner radius of the column and P(x)is the pressure at position xdescribed 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].


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 (nC14H30, nC16H32, nC20H42, nC25H52, nC30H62, nC35H72, nC40H82, nC50H102, nC60H122, nC70H142, and nC80H162).

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 “C15 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 C15 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): nC14, nC16, nC20, nC25, nC30, nC35, nC40, nC50, nC60, nC70, nC80.

  • 6 classes: alkene, CH4, C2H6, C3-C5, C6-C13 and C15 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 OilsnC14, nC16, nC20, nC25, nC30, nC35, nC40, nC45, nC50, nC55, nC60, nC65, nC70, nC75, nC80nC14, nC16, nC20, nC25, nC30, nC35, nC40, nC50, nC60, nC70, nC80

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 (nC12H26-nC98H196) [4].

  • the chemical reactions are described using kinetics and mechanistic modelling of 11 n-alkanes: nC14, nC16, nC20, nC25, nC30, nC35, nC40, nC50, nC60, nC70, nC80 and 6 class molecules: alkene, CH4, C2H6, C3-C5, C6-C13 and C15 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 nC25H52 to nC98H198 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 kis the retention factor.


Note here that kis 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. ro and wcorrespond to the inner radius of the column and the film thickness of the stationary phase Ds,and Dmcorrespond to the diffusion constant respectively in the stationary and mobile phase. vmcorresponds to the velocity of migration of the carrier gas.

The local zone variance (σx2) 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.


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.


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 nC12H26, predicted using an iterative analytic approach11 using MATLAB and solving the diffusive-convection equation by finite element in COMSOL. (Column dimensionsTable 3and temperature programmingTable 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 (xi+1), is calculated, using Snijders [23, 27] approach Eq. (8) (see ref. [19]), the distribution factor (K), and the phase ratio (β).


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 (nC14, nC16, nC20, nC25, nC30, nC35, nC40, nC50, nC60, nC70, nC80) 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 (nC14, nC16, nC20, nC25, nC30, nC35, nC40, nC50, nC60, nC70, nC80) 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 nC14H30 and nC16H34 with the temperature program used (Table 2), and their associated short residence times inside the GC column. Similarly, within the range nC20H42 to nC40H82, insignificant conversion occurs, whereas in the case of nC50H102 the maximum mass loss through thermal decomposition before elution is 0.003%.

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. nC60H122 has a significant loss in the stationary phase where only 2.43·10−12 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. nC70H142, presents a cumulative conversion of 0.001% at 385°C with only 2.32·10−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 nC60H122 decomposition starts, nC50H102 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 nC60H122 are released gradually, evidenced by a slowly increasing baseline signal.

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

Lastly, nC80H162 starts to decompose at 0.41 m from the column inlet, while nC70H142 is located 1.64 m from the inlet. Thus, when nC70H142 is essentially fully eluted (99.99%), at 7.83 m from the column inlet, nC80H162 has undergone a cumulative conversion of 0.52% mass loss by pyrolysis, relative to mass injected. That equates to 3.97.10−11 g of nC80H162 converted into pyrolysis products, and which co-elutes with nC70H142, 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 nC12H26 to nC98H198, [4] was used as main input for the calculation of the degree of elution of each of the n-alkanes studied.

The degree of elutionhas 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 C14H30 to nC80H162. Degree of elution = moles of “i” inside the GC column at time (t) /moles injected of “i”.

Alkanes heavier than nC60H122 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).

nC70H142 starts to elute at 29 minutes with a 99.99% degree of elution at 31.3 minutes, and 100% at 31.5 minutes. nC70H142 takes 2.5 minutes to elute completely.

nC80H162 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. nC80H162 takes 7.1 minutes to elute and 8.5 minutes to completely elute.

In this simulated study, components from nC70H142 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 nC80H162 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 (nC14H30, nC16H32, nC20H42, nC25H52, nC30H62, nC35H72, nC40H82, nC50H102, nC60H122, nC70H142, and nC80H162). 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 nC80H162, thermally decomposed before nC70H142. Therefore, co-elution of nC70H142 and the pyrolysis product of nC80H162 makes the GC analysis of nC70H142 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 nC70H142 take progressively longer to elute completely from the column, viz. nC70H142 takes 2.3 minutes and nC80H162 takes 7.1 minutes, with co-elution of decomposition products in each case compromising their analyses.

Finally, nC80H162 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 nC80H162 analysis may not be utilised under these HTGC conditions.



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.


DeffEffective average Diffusivity (unit length2/unit time)
DApparent diffusion coefficient, represents all factors causing dispersion (unit length2/unit time)
DMDiffusion constant, mobile phase (unit length2/unit time)
DSDiffusion constant, stationary phase (unit length2/unit time)
H(x,t))Column plate height, spatial rate of dispersion of a zone (unit length)
KDistribution factor of a compound (moles/volume) in stationary phase/(moles/volume) in gas phase)
kretention factor of a compound (moles in stationary phase/moles in gas phase).
LLength of the GC column (unit length)
m(x,to)mass profile for every analyte (particles/unit length).
Ni,MMoles of component “i” in the mobile phase.
Ni,SMoles of component “i” in the stationary phase.
P(x)Pressure at position x (Pa)
PinPressure at the GC colum inlet.
PoutPressure at the GC column outlet (Pa)
roInternal radius of GC column (unit length)
ramp TRamp of temperature of the temperature programmed.
T(t)Temperature at the time t.
ToInitial temperature of the temperature programmed.
ttime (unit time)
veffEffective cross-sectional average velocity (unit length/unit time)
vMVelocity of migration of the carrier gas (unit length/unit time)
wFilm thickness (unit length)
XiFraction of component i in the gas phase relative to the moles in both stationary and gas phase.
x0Centroid of Gaussian distribution of distribution of component inside the GC column (unit length)
xPosition 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).
ηmViscosity of the carrier gas.(μPa·s).
ΔtTime step (unit time).


  1. 1. MATLAB. version 7.10.0 (R2010a). Natick, Massachusetts: The MathWorks Inc.; 2010
  2. 2. Multiphysics C. Introduction to COMSOL multiphysics extregistered. COMSOL Multiphysics, Burlington, MA, accessed Feb. 1998;9:2018
  3. 3. CHEMKIN 10112, Reaction Design: San Diego, 2011
  4. 4. Hernandez-Baez DM, Reid A, Chapoy A, Tohidi B. Determination of distribution factors for heavy n-alkanes (nC12-nC98) in high temperature gas chromatography. J Chromatogr A. 2019;1591:138–146
  5. 5. SGE. HT5 GC Columns.
  6. 6. ASTM D7169-11. (Standard Test Method for boiling Point Distribution of Samples with Residues such as Crude Oils and Atmospheric and Vacuum Residues by High Temperature Gas Chromatography)
  7. 7. Hernandez-Baez DM, Reid A, Chapoy A, Tohidi B, Bounaceur R. Establishing the Maximum Carbon Number for Reliable Quantitative Gas Chromatographic Analysis of Heavy Ends Hydrocarbons. Part 3. Coupled Pyrolysis-GC Modeling. Energy Fuels. 2019 Mar 21;33(3):2045–2056
  8. 8. Keyes DE, Mcinnes LC, Woodward C, Gropp W, Myra E, Pernice M, Bell J, Brown J, Clo A, Connors J, Constantinescu E, Estep D, Evans K, Farhat C, Hakim A, Hammond G, Hansen G, Hill J, Isaac T, Jiao X, Jordan K, Kaushik D, Kaxiras E, Koniges A, Lee K, Lott A, Lu Q, Magerlein J, Maxwell R, Mccourt M, Mehl M, Pawlowski R, Randles AP, Reynolds D, Rivière B, Rüde U, Scheibe T, Shadid J, Sheehan B, Shephard M, Siegel A, Smith B, Tang X, Wilson C, Wohlmuth B. Multiphysics simulations. Int J High Perform Comput Appl. 2013;27(1):4–83
  9. 9. Giddings JC. Dynamics of chromatography, Par I, Principles and theory. Edw Arnold Publ Lond Marcel Dekker Inc N Y. 1965;
  10. 10. Gonzalez FR, Alessandrini JL, Nardillo AM. Revision of a theoretical expression for gas-liquid chromatographic retention. J Chromatogr A. 1999;852(2):583–588
  11. 11. Ben Naim A. Solvation Thermodyanmics. Plenum Press. 1987;New York
  12. 12. Gonzalez FR. Interpreting the gas chromatographic retention of n-alkanes. J Chromatogr A. 2000;873(2):209–219
  13. 13. Aldaeus F, Thewalim Y, Colmsjo A. Prediction of retention times and peak widths in temperature-programmed gas chromatography using the finite element method. J Chromatogr A. 2009 Jan;1216(1):134–139
  14. 14. Gonzalez FR. Considerations on the temperature dependence of the gas-liquid chromatographic retention. J Chromatogr A. 2002;942(1–2):211–221
  15. 15. Castells RC, Arancibia EL, Nardillo AM. Regression against temperature of gas-chromatographic retention data. J Chromatogr. 1990;504(1):45–53
  16. 16. Aldaeus F. New tools for trapping and separation in gas chromatgorahpy and dielectrophoresis (Improved performance by aid of computer simulation). Dr Thesis Anal Chem Stockh Univ. 2007;
  17. 17. Snijders H, Janssen HG, Cramers C. Optimization of temperature-programmed gas chromatographic separations .1. Prediction of retention times and peak widths from retention indices. J Chromatogr A. 1995;718(2):339–355
  18. 18. Aldaeus F, Thewalim Y, Colmsjo A. Prediction of retention times of polycyclic aromatic hydrocarbons and n-alkanes in temperature-programmed gas chromatography. Anal Bioanal Chem. 2007 Oct;389(3):941–950
  19. 19. Hernandez-Baez DM, Reid A, Chapoy A, Tohidi B, Bounaceur R. Establishing the Maximum Carbon Number for Reliable Quantitative Gas Chromatographic Analysis of Heavy Ends Hydrocarbons. Part 2. Migration and Separation Gas Chromatography Modeling. Energy Fuels. 2013;27(4):2336
  20. 20. Davankov VA. The true physical meaning of the corrected retention volumes in GC. Chromatographia. 1997;44(5–6):279–282
  21. 21. Kestin J, Knierim K, Mason EA, Najafi ST, Ro ST, Waldman M. Equilibrium and trasport properties of the noble gases and their mixtures at low density. J Phys Chem Ref Data. 1984;13(1):229
  22. 22. Hawkes SJ. Viscosities of carrier gases at gas-chromatograph temperatures and pressures. Chromatographia. 1993 Oct;37(7–8):399–401
  23. 23. Hernandez-Baez DM, Tohidi B, Chapoy A, Bounaceur R, Reid A. Establishing the Maximum Carbon Number for Reliable Quantitative Gas Chromatographic Analysis of Heavy Ends Hydrocarbons. Part 1: Low-Conversion Thermal Cracking Modeling. Energy Fuels. 2012 May 17;26(5):2600–2610
  24. 24. Golay MJE. Theory of Chromatography in Open and Coated Tubular Columns with Round and Rectangular Cross-Sections. 1958. 36 p
  25. 25. Blumberg LM. Temperature-Programmed Gas Chromatography. 2010
  26. 26. Giddings JC, Seager SL, Stucki LR, Stewart GH. Plate Height in gas chromatography. Anal Chem. 1960;32(8):867
  27. 27. Snijders H, Janssen HG, Cramers C. Optimization of temperature-programmed gas chromatographic separations 0.1. Prediction of retention times and peak widths from retention indices. J Chromatogr A. 1995;718(2):339

Written By

Diana Margarita Hernandez-Baez, Alastair Reid, Antonin Chapoy, Bahman Tohidi, Roda Bounaceur and François Montel

Submitted: January 22nd, 2021 Reviewed: May 28th, 2021 Published: August 5th, 2021