A Method for 3D Modeling of Chemical Separation in Microfabricated Gas Chromatography Columns: Case Studies in Temperature Gradients and Stationary Phase Topologies

: Conventional capillary gas chromatography (GC) columns, which have circular symmetry in cross-section and uniformity in length, are well modeled mathematically by the GC rate theory. However, even after adaptation, the theory has limited applicability to many unconventional properties in microfabricated GC columns, such as trapezoidal cross-sections, non-uniform stationary phase, and temperature gradients. This paper reports a 3D ﬁnite-element model for the chemical separation process in microfabricated GC columns using COMSOL. The model incorporates gas ﬂow, diffusion, partition, and temperature effects, enabling quantitative assessment of the separation performance of microfabricated GC columns with different stationary phase coating topologies and temperature gradients. To address the tremendous computational burden in such a 3D model, this paper investigates methods of providing proper meshing and dimensional scaling. For validation purposes, the implemented model was ﬁrst applied to a conventional capillary GC column and showed good matches to both the analytical calculation and experimental results. Next, the model was used to assess microfabricated columns with a trapezoidal cross-section and different stationary phase coating topologies. The results showed that, for the cases under consideration, a single-side-coated column provides only a 33% lower separation resolution compared to a double-side-coated column, and a parabolic stationary phase proﬁle provides only a 12% lower separation resolution compared to a uniform proﬁle. The model also indicated that temperature gradients have a negligible impact on separation performance.


Introduction
In analytical chemistry, gas chromatography (GC) remains a gold-standard technique for complex chemical analyses. GC separates complex mixtures into individual compounds, which can be detected using a detector or identified using a spectrometer. In GC, the separation column is a vital component. It is a gas flow channel coated with a stationary phase that provides separation of the compounds passing through the column. Depending on the chemical and physical properties, each compound establishes a certain level of partition between the stationary phase and the carrier gas, thereby eluting the column at a certain time, i.e., the retention time, which is characteristic of the compound [1]. GC has a multitude of applications, and more are emerging. Examples include the analysis of new products in the pharmaceuticals industry, the monitoring of metabolites in biological systems, environmental pollution monitoring, and the analysis of petroleum products and in the petroleum industry [1]. Although GC is powerful and versatile, conventional GC instruments are relatively bulky pieces of high-power benchtop equipment. Microscale gas chromatographs (µGCs) are miniaturized GC instruments with high portability and low power consumption, enabling timely on-site measurements. µGCs typically incorporate microfabricated components, such as microfabricated preconcentrators, columns, and detectors [2][3][4][5][6].
While the separation columns typically used in conventional GC instruments are capillary columns with circular cross-sections, the cross-section of a µGC column is determined by its microfabrication method and is often non-circular. For instance, a µGC column fabricated using deep reactive-ion etching (DRIE) typically has a rectangular cross-section. In some work, during DRIE, dense arrays of narrow and vertical pillars may also be created in the channel, forming a "semi-packed" column [7]. In comparison, a µGC column fabricated using wet etching, sandblasting, or metal micromachining techniques may have tapered sidewalls or curved surfaces [8,9]. Furthermore, the stationary phase topology in a column is determined by its coating process. Conventionally, separation columns are coated using static or dynamic methods, which typically produce a uniform coating on the entire inner wall of the capillary column. However, for microfabricated columns, these methods typically leave a thicker coating at the corners than on the walls, i.e., the pooling effect, which causes undesired peak broadening [6]. Additionally, these methods have limited throughputs. Due to these challenges, alternative coating methods are being explored, such as chemical vapor deposition [10,11], atomic layer deposition [12][13][14], and spin-coating [15]. These alternative methods may only coat a portion of the inner surface of the column.
When designing and developing new fabrication and coating methods for µGC columns, it is highly valuable to obtain a theoretical prediction of the separation performance before delving into the labor-intensive and time-consuming experimental practice. The theory is well established for columns with circular or rectangular cross-sections and a stationary phase uniformly coated on the entire inner wall [16,17]. Mathematical formulae for columns with rectangular cross-sections and one or two walls coated are also reported [18]. However, theoretical models are not yet available for other cross-sections (e.g., semi-circular or semi-trapezoidal) and other stationary phase topologies (e.g., non-uniform or patterned coating on selected walls), which require 3D modeling.
Although the impact of column temperature has been well studied, the established theory mostly assumes a uniform column temperature that is either constantly maintained or programmed. However, the impact of a temperature gradient within the separation column is also worth investigating. Historically, it has been under debate whether separation performance can be improved by a longitudinal temperature gradient, i.e., with the upstream end of the column at a higher temperature than the downstream end. Some have argued that, when flowing along such a temperature gradient, a peak would tend to become sharper, as the trailing portion of the peak would be at a higher temperature than the leading portion and, therefore, would tend to catch up with the latter [19][20][21][22][23]. Others have refuted that such a temperature gradient would not enhance separation, as the peaks behind would have the same tendency to catch up with the peaks ahead [24][25][26][27][28]. Additionally, within a practical µGC system, the temperature distribution of the separation column may be affected by the package and thermal crosstalk from other components, forming a temperature gradient that possibly does not align with the flow direction. A model that can incorporate temperature effects can provide insight into these problems.
This work aims to develop a methodology to model the µGC column using COMSOL. It is intended that the developed model will enable the quantitative evaluation of different stationary phase topologies. It is also intended that temperature effects will be incorporated, enabling the investigation of temperature gradients on separation performance. The model can also allow visualization of the chemical concentration profile travelling along the column, providing an exclusive insight into the separation process.
To meet the goals of this work, a 3D FEA model is needed, which poses challenges associated with the substantial differences in various dimensions of the column. For a good mesh in the FEA, the aspect ratio (i.e., the ratio of the largest dimension to the smallest dimension) of the meshed elements is preferably close to unity. However, in a µGC column, the length is often 3-4 orders of magnitude larger than the cross-sectional dimensions, which are often another 2-3 orders of magnitude larger than the stationary phase thickness. Consequently, a tremendously large number of meshed elements is needed to provide the desired accuracy of the result, which is time consuming. To address this challenge, an artificially thick stationary phase is set with scaled material properties to represent a thinner layer, thereby significantly reducing the required number of mesh elements and the simulation time. The methodology to model the separation column is described in Section 2 and the results and discussion are presented in Section 3, followed by the conclusions in Section 4. For all the figures shown in this paper, the raw images were exported from COMSOL ® and MATLAB ® , and subsequently labeled in Inkscape ® 1.1.

Modeled Geometries
The objectives of the model include the study of the microfabricated columns with a trapezoidal cross-section, and further, the study of the stationary phase topologies therein. To validate the modeling approach, a straight capillary column is first modeled and compared to the experimental results. The modeled capillary column has a length of 1 m, an inner diameter of 250 µm, and a uniformly coated polydimethylsiloxane (PDMS) stationary phase with a thickness of 1 µm (Figure 1a).
provide the desired accuracy of the result, which is time consuming. To address this challenge, an artificially thick stationary phase is set with scaled material properties to represent a thinner layer, thereby significantly reducing the required number of mesh elements and the simulation time. The methodology to model the separation column is described in Section 2 and the results and discussion are presented in Section 3, followed by the conclusions in Section 4. For all the figures shown in this paper, the raw images were exported from COMSOL ® and MATLAB ® , and subsequently labeled in Inkscape ® 1.1.

Modeled Geometries
The objectives of the model include the study of the microfabricated columns with a trapezoidal cross-section, and further, the study of the stationary phase topologies therein. To validate the modeling approach, a straight capillary column is first modeled and compared to the experimental results. The modeled capillary column has a length of 1 m, an inner diameter of 250 μm, and a uniformly coated polydimethylsiloxane (PDMS) stationary phase with a thickness of 1 μm (Figure 1a).
The modeled microfabricated column consists of 38 straight segments connected via semi-circular segments (Figure 1b). Each straight segment is 13.6 mm long. The cross-section of the microfabricated column is trapezoidal with a 500 μm-wide bottom edge, a 371 μm-wide top edge, a 160 μm height, and two sidewalls with a 22° inward taper. For the microfabricated column, three different cases of the stationary phase topology are modeled. In Case 1, the stationary phase has parabolic profiles on both the top and bottom surfaces of the column, with a 5 μm thickness at the center ( Figure 1c). The widths of the stationary phase on the top and bottom surfaces are 371 μm and 400 μm, respectively. In Case 2, the stationary phase is uniformly coated on both the top and bottom surfaces of the column and has a thickness of 2.95 μm (Figure 1d). In Case 3, the stationary phase is uniformly coated only on the bottom surface of the column and has a thickness of 5.14 μm (Figure 1e). All three cases have equal volumes of the stationary phase.  The modeled microfabricated column consists of 38 straight segments connected via semi-circular segments (Figure 1b). Each straight segment is 13.6 mm long. The crosssection of the microfabricated column is trapezoidal with a 500 µm-wide bottom edge, a 371 µm-wide top edge, a 160 µm height, and two sidewalls with a 22 • inward taper. For the microfabricated column, three different cases of the stationary phase topology are modeled. In Case 1, the stationary phase has parabolic profiles on both the top and bottom surfaces of the column, with a 5 µm thickness at the center (Figure 1c). The widths of the stationary phase on the top and bottom surfaces are 371 µm and 400 µm, respectively. In Case 2, the stationary phase is uniformly coated on both the top and bottom surfaces of the column and has a thickness of 2.95 µm (Figure 1d). In Case 3, the stationary phase is uniformly coated only on the bottom surface of the column and has a thickness of 5.14 µm (Figure 1e). All three cases have equal volumes of the stationary phase.
The selected thicknesses of the stationary phase noted above are modeled using a hypothetical thickness of 20 µm COMSOL geometry. This hypothetical thickness reduces the number of elements in the mesh and avoids re-meshing the model each time a different thickness of the stationary phase needs to be studied. A scaling factor n (explained in Section 2.3) is employed in the simulation to obtain the results for the actual thickness.

Modeling of the Gas Flow
The gas flow regime in a separation column is determined using the Reynolds number (Re), which can be calculated as: where ρ and µ are the density and dynamic viscosity of the carrier gas, respectively; U is the characteristic velocity, i.e., the average gas flow speed; and d H is the characteristic length, i.e., the hydraulic diameter of the column cross-section [29]. For air at standard temperature and pressure, ρ = 1. Therefore, this work implements the "Laminar Flow" module in COMSOL to model the incompressible flow of carrier gas in the separation column, which applies the following governing equations: where u is the velocity field of the gas, p is the pressure, = τ is the stress tensor, and → f is the body force.
For the 'Laminar Flow' module, the model domain is initialized with zero velocity, and a no-slip boundary condition is applied to the column walls, specifying that the gas flow in contact with the column walls has zero velocity. The inlet boundary condition is set with a fully developed flow at an average velocity U (with the value indicated above). The outlet boundary condition is set at 0 Pa relative to a reference pressure at 1 atm.

Modeling of the Chemical Transport
The chemicals transported along the separation column undergo convection, diffusion, and partition between the carrier gas and the stationary phase. These transportation mechanisms are modeled by the COMSOL module "Transport of Diluted Species" in this work.
Equations (5)-(7) model the convection (coupled with the flow field) and diffusion in the carrier gas, the diffusion in the stationary phase, and the partition between the two phases, respectively. Here, c is the concentration of the chemical, t is time, D is the diffusivity of the chemical, S is a source/sink term, and K is the partition coefficient. The subscripts CG and SP denote the carrier gas and the stationary phase, respectively.
Because this work scales the stationary phase thickness to a large, fixed hypothetical value (20 µm) to address the meshing and geometrical challenges, a scaling factor n, as mentioned in Section 2.1, must be applied to compensate for the impacts of such scaling. The n is defined as the ratio of actual stationary phase thickness to the hypothetical 20 µm. As a result, diffusivity in the stationary phase and the partition coefficient must be scaled as follows: Note that such scaling is applied only in the direction of the stationary phase thickness, i.e., only in the xand y-directions for the modeled capillary column and only in the ydirection for the modeled microfabricated columns. In COMSOL, such directional scaling is implemented by setting the diffusivity in the stationary phase as an anisotropic tensor.
For the 'Transport of Diluted Species' module, a no-flux boundary condition is applied to the column walls, specifying that no chemicals can penetrate through the column walls. The partition equilibrium condition is also applied to the interface between the stationary phase and carrier gas ( Table 1). The model domain is initialized with zero chemical concentration. To model the chemical injection at the inlet, a Gaussian concentration profile of the chemical is set to pass the inlet cross-section at 1 s, with a peak width at half-height (PWHH inlet ) of 0.5 s and a peak concentration of 4.4 × 10 −5 mol/m 3 .

Modeling of the Temperature Effects
The impacts of temperature on the separation result from the temperature dependence of the diffusivity and the partition coefficient. For a chemical α in another chemical β, the diffusivity (D α/β ) is a function of temperature T (in the unit of K) [30]: where M denotes the molecular weight of the chemical, p is pressure in atm, and ∑v i denotes the diffusion volume of the chemical. This work selects air as the carrier gas, which has a diffusion volume of 20.1 cm 3 [30]. This work selects hexane and benzene as the two representative chemicals for the model. Based on the known values of D hexane/air and D benzene/air at 25 • C [31], the diffusion volumes of hexane and benzene are extracted to be 140.88 cm 3 and 82.16 cm 3 , respectively.
In the stationary phase, the temperature dependence of the diffusivity can generally be described by the Arrhenius equation [32], D α/β = D 0 exp(−E A /RT), where R is the universal gas constant, D 0 is the maximal diffusion coefficient at infinite temperature, and E A is the activation energy for diffusion. However, the values of D 0 and E A are not always available. Therefore, in this work, for simplicity, the diffusivity values of hexane and benzene in PDMS (D hexane/PDMS and D benzene/PDMS , respectively) are approximated as temperature-independent constants at 2.5 × 10 −10 m 2 /s [33] and 2.8 × 10 −10 m 2 /s [34], respectively. This approximation holds well for the small temperature range of 20-40 • C investigated in this work.
The partition coefficient K as a function of temperature (in K) can be calculated from its known value K T 0 at a reference temperature T 0 (which is usually 298.15 K) [35].
where A, B, and C are Antoine constants (Table 1).

Simulation Setup
In the COMSOL simulation, two Study steps are sequentially performed. The first step models the gas flow; this is performed as a Stationary Study, because the flow is fully developed and remains constant during the separation. The second step models the chemical transport process; this is performed as a Time-Dependent Study, because the chemical concentrations continue to change during the separation.
In the second step, because the PWHH inlet is set to 0.5 s, the time-step used during 0-2 s is set to 0.01 s to provide a smooth inlet Gaussian profile. Beyond 2 s, a time-step of 0.02 s is selected to reduce the simulation time by half without significantly impacting the result. These values are obtained from the mesh convergence study (Section 3.3).
The impact of the longitudinal temperature gradient is investigated using the modeled capillary column. Three temperature cases are considered: along the column length (in the z-direction) from the inlet to the outlet, Case I applies a temperature distribution linearly decreasing from 40 • C to 20 • C, Case II applies a temperature distribution linearly increasing from 20 • C to 40 • C, and Case III applies a uniform 30 • C ( Figure 2a).
The partition coefficient K as a function of temperature (in K) can be calculated from its known value at a reference temperature T0 (which is usually 298.15 K) [35].
where A, B, and C are Antoine constants (Table 1).

Simulation Setup
In the COMSOL simulation, two Study steps are sequentially performed. The first step models the gas flow; this is performed as a Stationary Study, because the flow is fully developed and remains constant during the separation. The second step models the chemical transport process; this is performed as a Time-Dependent Study, because the chemical concentrations continue to change during the separation.
In the second step, because the PWHHinlet is set to 0.5 s, the time-step used during 0-2 s is set to 0.01 s to provide a smooth inlet Gaussian profile. Beyond 2 s, a time-step of 0.02 s is selected to reduce the simulation time by half without significantly impacting the result. These values are obtained from the mesh convergence study (Section 3.3).
The impact of the longitudinal temperature gradient is investigated using the modeled capillary column. Three temperature cases are considered: along the column length (in the z-direction) from the inlet to the outlet, Case I applies a temperature distribution linearly decreasing from 40 °C to 20 °C, Case II applies a temperature distribution linearly increasing from 20 °C to 40 °C, and Case III applies a uniform 30 °C (Figure 2a).
The impact of the lateral temperature gradient is investigated using the modeled microfabricated column. Similar to the capillary column model, three temperature cases are also considered: along the column width (in the x-direction) from the inlet to the outlet, Case IV applies a temperature distribution linearly decreasing from 40 °C to 20 °C, Case V applies a temperature distribution linearly increasing from 20 °C to 40 °C, and Case VI applies a uniform 30 °C (Figure 2b).  The impact of the lateral temperature gradient is investigated using the modeled microfabricated column. Similar to the capillary column model, three temperature cases are also considered: along the column width (in the x-direction) from the inlet to the outlet, Case IV applies a temperature distribution linearly decreasing from 40 • C to 20 • C, Case V applies a temperature distribution linearly increasing from 20 • C to 40 • C, and Case VI applies a uniform 30 • C (Figure 2b).

Mesh Development
For the capillary column, a user-controlled mesh is applied to both study steps. In this approach, the cross-section of the column is first meshed with 2D elements (Figure 3a). This mesh is then swept along the length of the column to mesh the entire column geometry (Figure 3b). Such a mesh conforms well to the circular cross-section and enables a mesh convergence study of the impact of mesh element size along the column length.

Mesh Development
For the capillary column, a user-controlled mesh is applied to both study steps. In this approach, the cross-section of the column is first meshed with 2D elements ( Figure  3a). This mesh is then swept along the length of the column to mesh the entire column geometry (Figure 3b). Such a mesh conforms well to the circular cross-section and enables a mesh convergence study of the impact of mesh element size along the column length.
For the microfabricated column, a physics-controlled mesh, which is found to represent the column dimensions well, is used for the Stationary Study of gas flow. Next, a user-controlled mesh is used for the Time-Dependent Study, similar to the capillary column mesh (Figure 4).

Results and Discussion
This section describes the results obtained from the models. First, the gas flow velocity profiles and chemical concentration profiles during transport are presented. Then, the results of the mesh convergence study are presented and compared to the analytical calculation and experimental results. Next, the performance of the modeled microfabricated columns with different stationary phase topologies are compared. Finally, the impacts of the different temperature gradients on the separation columns are evaluated. Figure 5 shows the velocity distribution of the carrier gas in both the modeled capillary column and microfabricated column. As expected, the velocity in each column reached its maximum at the center of the cross-section, and was reduced to zero at the For the microfabricated column, a physics-controlled mesh, which is found to represent the column dimensions well, is used for the Stationary Study of gas flow. Next, a usercontrolled mesh is used for the Time-Dependent Study, similar to the capillary column mesh (Figure 4).

Mesh Development
For the capillary column, a user-controlled mesh is applied to both study steps. In this approach, the cross-section of the column is first meshed with 2D elements ( Figure  3a). This mesh is then swept along the length of the column to mesh the entire column geometry (Figure 3b). Such a mesh conforms well to the circular cross-section and enables a mesh convergence study of the impact of mesh element size along the column length.
For the microfabricated column, a physics-controlled mesh, which is found to represent the column dimensions well, is used for the Stationary Study of gas flow. Next, a user-controlled mesh is used for the Time-Dependent Study, similar to the capillary column mesh (Figure 4).

Results and Discussion
This section describes the results obtained from the models. First, the gas flow velocity profiles and chemical concentration profiles during transport are presented. Then, the results of the mesh convergence study are presented and compared to the analytical calculation and experimental results. Next, the performance of the modeled microfabricated columns with different stationary phase topologies are compared. Finally, the impacts of the different temperature gradients on the separation columns are evaluated. Figure 5 shows the velocity distribution of the carrier gas in both the modeled capillary column and microfabricated column. As expected, the velocity in each column reached its maximum at the center of the cross-section, and was reduced to zero at the

Results and Discussion
This section describes the results obtained from the models. First, the gas flow velocity profiles and chemical concentration profiles during transport are presented. Then, the results of the mesh convergence study are presented and compared to the analytical calculation and experimental results. Next, the performance of the modeled microfabricated columns with different stationary phase topologies are compared. Finally, the impacts of the different temperature gradients on the separation columns are evaluated. Figure 5 shows the velocity distribution of the carrier gas in both the modeled capillary column and microfabricated column. As expected, the velocity in each column reached its

Concentration Profile
The chemical separation process can be visualized from the results. For example, Figure  6 illustrates the separation process of hexane and benzene in the modeled microfabricated column of Case 3 operated at 25 °C. After the injection from the column inlet in the lower left corner, at t = 2 s, the profiles of hexane and benzene remained largely overlapping ( Figure 6a). As these profiles moved along the column, they exhibited partial separation at 10 s ( Figure  6b), well-established separation at 20 s (Figure 6c), followed by the full elution of hexane by 30 s, upon which only benzene remained in the column (Figure 6d).

Mesh Convergence and Experimental Verification
To ensure the accuracy of modeling, the mesh must be fine enough that the results converge, upon which a finer mesh does not significantly alter the results. Because the column length is at least 3-4 orders of magnitude larger than any dimension on the crosssection, it is critical to evaluate the impact of mesh fineness along the column length. Therefore, in the mesh convergence study performed in this work, a fixed mesh of 162 elements was applied to the cross-section of the column, whereas the mesh element count per meter (nL) of the column length was evaluated over a range of 2000-8000. Additionally, the impact of the time-step on the convergence was evaluated. This was performed by

Concentration Profile
The chemical separation process can be visualized from the results. For example, Figure 6 illustrates the separation process of hexane and benzene in the modeled microfabricated column of Case 3 operated at 25 • C. After the injection from the column inlet in the lower left corner, at t = 2 s, the profiles of hexane and benzene remained largely overlapping ( Figure 6a). As these profiles moved along the column, they exhibited partial separation at 10 s (Figure 6b

Concentration Profile
The chemical separation process can be visualized from the results. For example, Figure  6 illustrates the separation process of hexane and benzene in the modeled microfabricated column of Case 3 operated at 25 °C. After the injection from the column inlet in the lower left corner, at t = 2 s, the profiles of hexane and benzene remained largely overlapping ( Figure 6a). As these profiles moved along the column, they exhibited partial separation at 10 s ( Figure  6b), well-established separation at 20 s (Figure 6c), followed by the full elution of hexane by 30 s, upon which only benzene remained in the column (Figure 6d).

Mesh Convergence and Experimental Verification
To ensure the accuracy of modeling, the mesh must be fine enough that the results converge, upon which a finer mesh does not significantly alter the results. Because the column length is at least 3-4 orders of magnitude larger than any dimension on the crosssection, it is critical to evaluate the impact of mesh fineness along the column length. Therefore, in the mesh convergence study performed in this work, a fixed mesh of 162 elements was applied to the cross-section of the column, whereas the mesh element count per meter (nL) of the column length was evaluated over a range of 2000-8000. Additionally, the impact of the time-step on the convergence was evaluated. This was performed by

Mesh Convergence and Experimental Verification
To ensure the accuracy of modeling, the mesh must be fine enough that the results converge, upon which a finer mesh does not significantly alter the results. Because the column length is at least 3-4 orders of magnitude larger than any dimension on the cross-section, it is critical to evaluate the impact of mesh fineness along the column length. Therefore, in the mesh convergence study performed in this work, a fixed mesh of 162 elements was applied to the cross-section of the column, whereas the mesh element count per meter (n L ) of the column length was evaluated over a range of 2000-8000. Additionally, the impact of the time-step on the convergence was evaluated. This was performed by varying the time-step at certain fixed n L values. This section illustrates the mesh convergence study results for the capillary column at 25 • C, and a similar study can be also performed on the microfabricated column.
To evaluate the mesh convergence, four performance parameters for the column were extracted from the results and considered ( Table 2). These include (1) the retention time (t R ), which was relative to the injection time at t = 1 s; (2) the peak width at half-height (PWHH) of the chemical peaks eluting from the column; (3) the number of theoretical plates (N), which is an indicator of the column efficiency and is computed as t R PW HH 2 (12) and (4) the separation resolution (R s ), which is a measure of how well two peaks are separated and is computed as [16] R S = 1.18 t R.benzene − t R.hexane PW HH benzene + PW HH hexane (13)  While all four performance parameters showed trends of convergence, the PWHH benzene in Table 2 was selected as an example and plotted against n L (Figure 7). As is evident from Figure 7, the PWHH benzene exhibited an exponential trend with respect to n L and a projected converged value of 1.33 s. For each simulated Mesh Case, the error of the results in % was calculated with respect to the converged value ( Table 2). The case with the finest mesh, i.e., Mesh Case xv, had an error of only 6.77%, which was considered accurate enough for this work. The model was verified using both an analytical calculation approach and an experimental approach, both applied to the modeled capillary column. The analytical calculation approach was based on the standard GC plate theory and rate theory, which are well established for capillary columns. In this approach, the tR, height equivalent to a theoretical plate (HETP), and N can be estimated as follows [1,16] = (1 + ) = 2 + 1 + 6 + 11 96(1 + ) where k is the retention factor, which is calculated as k = K/βP, in which βP is the phase ratio, i.e., the volume ratio between the carrier gas and the stationary phase; L, D, and h are the column length, inner diameter, and stationary phase thickness, respectively; DCG and DSP are the diffusivities of the chemical in the carrier gas and in the stationary phase, respectively; is the variance of a Gaussian peak introduced by the extra-column effects, which represents the injection peak width in the model. Numerically, for a Gaussian injection peak, = /2.355. After obtaining the results from Equations (14) and (16), the PWHHhexane, PWHHbenzene, and RS can be computed from Equations (12) and (13).
The results of this analytical calculation approach are listed in the second-to-last row in Table 2. Compared to the analytical calculation results, Mesh Case xv showed <22% differences in tR and PWHH, and <13% differences in N and RS.
In the experimental approach to verifying the model, a commercial off-the-shelf (COTS) capillary column (Part # 115-2504, Ohio Valley Specialty Company, Marietta, OH, USA) of equivalent properties to the model was experimentally tested in a benchtop GC system (Model # 6890, Agilent Technologies, Inc., Santa Clara, CA, USA), which provided an equivalent temperature and flow rate to the model. Samples of hexane and benzene were injected at the GC sample inlet and detected using a flame ionization detector (FID) after eluting the column. The experimental results were compared to the modeling result obtained from Mesh Case xv (Table 2 and Figure 8). The comparison showed tR differences at 11.03% and 1.68% for hexane and benzene, respectively, and PWHH differences at 28.79% and 8.97% for hexane and benzene, respectively. Note that for Figure 7, Mesh Cases (iii), (vii), and (xii) were excluded. These three cases were inserted to evaluate the impact of the time-step at a given n L . Based on the comparison between Mesh Cases iii vs. iv; viii vs. viii; and xii vs. xiii, the selected time-steps in these cases were shown to be fine enough without significantly affecting the results. For example, compared to Mesh Case xii, Mesh Case xiii reduced the time-steps by half, but only changed the PWHH benzene by 1.3%.
The model was verified using both an analytical calculation approach and an experimental approach, both applied to the modeled capillary column. The analytical calculation approach was based on the standard GC plate theory and rate theory, which are well established for capillary columns. In this approach, the t R , height equivalent to a theoretical plate (HETP), and N can be estimated as follows [1,16] where k is the retention factor, which is calculated as k = K/β P , in which β P is the phase ratio, i.e., the volume ratio between the carrier gas and the stationary phase; L, D, and h are the column length, inner diameter, and stationary phase thickness, respectively; D CG and D SP are the diffusivities of the chemical in the carrier gas and in the stationary phase, respectively; σ 2 ec is the variance of a Gaussian peak introduced by the extra-column effects, which represents the injection peak width in the model. Numerically, for a Gaussian injection peak, σ ec = PW HH inlet /2.355. After obtaining the results from Equations (14) and (16), the PWHH hexane , PWHH benzene , and R S can be computed from Equations (12) and (13).
The results of this analytical calculation approach are listed in the second-to-last row in Table 2. Compared to the analytical calculation results, Mesh Case xv showed <22% differences in t R and PWHH, and <13% differences in N and R S .
In the experimental approach to verifying the model, a commercial off-the-shelf (COTS) capillary column (Part # 115-2504, Ohio Valley Specialty Company, Marietta, OH, USA) of equivalent properties to the model was experimentally tested in a benchtop GC system (Model # 6890, Agilent Technologies, Inc., Santa Clara, CA, USA), which provided an equivalent temperature and flow rate to the model. Samples of hexane and benzene were injected at the GC sample inlet and detected using a flame ionization detector (FID) after eluting the column. The experimental results were compared to the modeling result obtained from Mesh Case xv (Table 2 and Figure 8). The comparison showed t R differences at 11.03% and 1.68% for hexane and benzene, respectively, and PWHH differences at 28.79% and 8.97% for hexane and benzene, respectively. Note that all the parameter values used in the model were obtained from the literature without being intentionally adjusted to match the analytical calculation or experimental results. The modest differences between the modeling and analytical calculation results are not atypical of other comparisons between analytical and modeling results; for example, analytical results generally assume idealized geometries and boundary conditions that are not discretized. The differences between the modeling and experimental results are expected from factors such as the manufacturing variations in the COTS capillary column, inaccuracies in the experimental flow rate and temperature, and minute variations in the chemical injection. Given the overall match to the analytical calculation and experimental results, the model is sufficiently accurate to investigate the impact of the stationary phase topologies and temperature gradients.

Stationary Phase Topologies in Microfabricated Columns
Based on the mesh convergence study described in Section 3.3, the three PDMS topologies cases described in Figure 1 were modeled and compared. A converged mesh with an nL of approximately 8800 (120 longitudinal elements in each straight segment and 20 elements in each curved segment of the microfabricated column) was applied to all three cases. In the modeling results, all three cases provided similar tR values, but noticeably different PWHH, N, and RS values (Figure 9 and Table 3). Compared to Case 2, Case 1 provided 19% broader hexane and benzene peaks, whereas Case 3 provided a 55% broader hexane peak and a 48% broader benzene peak. Consequently, for the separation performance, as represented by the resolution between hexane and benzene (denoted by RS), Case 1 and Case 3 were 12% and 33% worse than Case 2, respectively.
The similarity in the tR values was expected, as the overall volumes of PDMS in the three cases were the same. The overall worse performance of Case 3 than Case 2 was also expected because of two factors. First, the PDMS in Case 3 was only on the bottom surface, causing a larger effective distance for chemical diffusion from the carrier gas to its interface with the PDMS, i.e., causing larger resistance to mass transfer in the carrier gas [1]. Second, the PDMS in Case 3 was thicker, causing a larger effective distance for chemical diffusion from the PDMS to its interface with the carrier gas, i.e., causing larger resistance to mass transfer in the PDMS [1]. Nevertheless, compared to Case 2, the separation performance reduction in Case 3 may be considered tolerable for certain applications, especially when the benefit of simpler fabrication provided by Case 3 is emphasized.
Compared to Case 2, the separation performance reduction in Case 1 appears insignificant. This may be an encouraging sign for certain stationary phase coating processes, Note that all the parameter values used in the model were obtained from the literature without being intentionally adjusted to match the analytical calculation or experimental results. The modest differences between the modeling and analytical calculation results are not atypical of other comparisons between analytical and modeling results; for example, analytical results generally assume idealized geometries and boundary conditions that are not discretized. The differences between the modeling and experimental results are expected from factors such as the manufacturing variations in the COTS capillary column, inaccuracies in the experimental flow rate and temperature, and minute variations in the chemical injection. Given the overall match to the analytical calculation and experimental results, the model is sufficiently accurate to investigate the impact of the stationary phase topologies and temperature gradients.

Stationary Phase Topologies in Microfabricated Columns
Based on the mesh convergence study described in Section 3.3, the three PDMS topologies cases described in Figure 1 were modeled and compared. A converged mesh with an n L of approximately 8800 (120 longitudinal elements in each straight segment and 20 elements in each curved segment of the microfabricated column) was applied to all three cases. In the modeling results, all three cases provided similar t R values, but noticeably different PWHH, N, and R S values ( Figure 9 and Table 3). Compared to Case 2, Case 1 provided 19% broader hexane and benzene peaks, whereas Case 3 provided a 55% broader hexane peak and a 48% broader benzene peak. Consequently, for the separation performance, as represented by the resolution between hexane and benzene (denoted by R S ), Case 1 and Case 3 were 12% and 33% worse than Case 2, respectively.
The similarity in the t R values was expected, as the overall volumes of PDMS in the three cases were the same. The overall worse performance of Case 3 than Case 2 was also expected because of two factors. First, the PDMS in Case 3 was only on the bottom surface, causing a larger effective distance for chemical diffusion from the carrier gas to its interface with the PDMS, i.e., causing larger resistance to mass transfer in the carrier gas [1]. Second, the PDMS in Case 3 was thicker, causing a larger effective distance for chemical diffusion from the PDMS to its interface with the carrier gas, i.e., causing larger resistance to mass transfer in the PDMS [1]. Nevertheless, compared to Case 2, the separation performance reduction in Case 3 may be considered tolerable for certain applications, especially when the benefit of simpler fabrication provided by Case 3 is emphasized.
Micro 2022, 2, FOR PEER REVIEW 12 where a coated and patterned stationary phase film may reflow into a parabolic profile before being fully crosslinked and immobilized.

Impact of Temperature Gradients
The three cases of longitudinal temperature gradients in a capillary column (as described in Section 2.5) were modeled. The resulting chromatograms are overlaid in Figure  10, and the performance metrics are compared in Table 4. From these results, it is apparent that the negative longitudinal temperature gradient (in Case I) and the positive one (in Case II) provided equal performance, which was almost equal to the isothermal Case III. This result opposes the idea that the negative thermal gradient would significantly improve separation [19][20][21][22][23] and corroborates the counter-argument [24][25][26][27][28].   Compared to Case 2, the separation performance reduction in Case 1 appears insignificant. This may be an encouraging sign for certain stationary phase coating processes, where a coated and patterned stationary phase film may reflow into a parabolic profile before being fully crosslinked and immobilized.

Impact of Temperature Gradients
The three cases of longitudinal temperature gradients in a capillary column (as described in Section 2.5) were modeled. The resulting chromatograms are overlaid in Figure 10, and the performance metrics are compared in Table 4. From these results, it is apparent that the negative longitudinal temperature gradient (in Case I) and the positive one (in Case II) provided equal performance, which was almost equal to the isothermal Case III. This result opposes the idea that the negative thermal gradient would significantly improve separation [19][20][21][22][23] and corroborates the counter-argument [24][25][26][27][28].
For the three lateral temperature gradient cases in the modeled microfabricated column as described in Section 2.5, the applied temperature gradient was effectively perpendicular to the dominating flow direction. From the modeling results ( Figure 11 and Table 5), the negative lateral temperature gradient (in Case IV), the positive one (in Case V), and the isothermal Case VI appeared to provide almost the same performance (with <2% variation).
Such a result provides an encouraging and quantitative understanding of tolerance to temperature non-uniformity in the microfabricated columns. For the three lateral temperature gradient cases in the modeled microfabricated column as described in Section 2.5, the applied temperature gradient was effectively perpendicular to the dominating flow direction. From the modeling results ( Figure 11 and Table  5), the negative lateral temperature gradient (in Case IV), the positive one (in Case V), and the isothermal Case VI appeared to provide almost the same performance (with <2% variation). Such a result provides an encouraging and quantitative understanding of tolerance to temperature non-uniformity in the microfabricated columns.    For the three lateral temperature gradient cases in the modeled microfabricated column as described in Section 2.5, the applied temperature gradient was effectively perpendicular to the dominating flow direction. From the modeling results ( Figure 11 and Table  5), the negative lateral temperature gradient (in Case IV), the positive one (in Case V), and the isothermal Case VI appeared to provide almost the same performance (with <2% variation). Such a result provides an encouraging and quantitative understanding of tolerance to temperature non-uniformity in the microfabricated columns.    It is also worth noting that the isothermal cases always showed slightly smaller t R values than the corresponding cases with temperature gradients (Tables 4 and 5 and Figures 10 and 11). This is mathematically reasonable. As described by Equation (11), K(T) exhibits a convex function of T in the modeled temperature range of 20-40 • C. According to Jensen's inequality [38], K (30 • C) is smaller than the average of K(T) over 20-40 • C. Additionally, for a given column geometry, t R increases linearly with K, as indicated by Equation (14). Therefore, the t R for isothermal 30 • C is expected to be smaller than that with a constant temperature gradient over 20-40 • C.

Conclusions
This work presents a methodology for modeling µGC separation columns using 3D FEA, which is essential for columns with unconventional and asymmetric cross-sections and stationary phase coating topologies. The methodology starts with essential simplifications of the governing equations. Meshing of the separation column is performed in a customdefined manner to reduce simulation time while preserving enough accuracy of the results. A scaling factor of the stationary phase thickness is implemented to overcome the challenges in the large mesh density caused by the large span of column dimensions. Additionally, this scaling factor eliminates the need to re-mesh when modeling different thicknesses of the stationary phase with a given topology. A mesh convergence study is performed to confirm the required mesh density. The methodology is validated by the overall match between the modeling results of a capillary column and both the analytical calculation and experimental results.
Using the developed model, the impacts of stationary phase topology and temperature gradients on the performance of the separation column are quantitatively investigated. While a stationary phase coating on fewer inner surfaces or with a reflowed profile causes a reduction in separation performance, such a limited reduction may be well accepted if the benefits posed by the simpler fabrication process of the column are more desired. Compared to the isothermal cases, neither longitudinal nor lateral temperature gradients within the ranges considered in this work appear to affect the separation performance.
The developed model can be readily adapted to model various cross-sections of the column and topologies of its stationary phase. For example, it can be adapted to evaluate the impact of the pooling effect, which is common in microfabricated columns that use static or dynamic coating methods. Additionally, the model can assist in the exploration of any attractive stationary phase topologies. In the future, the model can be expanded to examine the impact of other µGC components on the overall separation performance of the system, e.g., the impact of a broad injection peak from the preconcentrator, particularly if the column is short; the impact of the fluidic pathway of the detector; and the impact of dead volume associated with other fluidic connections at the pump, valves, etc. The model can also be enhanced to incorporate adsorption at the column walls for certain surface-adsorptive analyte chemicals, such as amines and phosphonate esters.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.