Development and Validation of Retention Models in Supercritical Fluid Chromatography for Impregnation Process Design

Abstract

The retention factor is the key quantity for the thermodynamic analysis of the retention mechanism in chromatographic experiments. In this work, we measure retention factors for moderately polar solutes on four silica-based porous matrices as stationary phases by supercritical fluid chromatography. Elution of the solutes is only possible with binary mixtures of supercritical carbon dioxide (sc-CO₂) and modifier (methanol) due to the low polarity of pure sc-CO₂. The addition of modifiers makes the retention mechanism more complex and masks interactions between solute and stationary phase. In this work, we develop and validate several retention models that allow the obtaining of retention factors in modifier-free sc-CO₂. Such models pave the way for quantifying adsorption interactions between polar solutes and non-swellable porous matrices in pure sc-CO₂ based on retention data obtained in sc-CO₂/modifier mixtures. The obtained information will thereby facilitate the understanding and design of impregnation processes, which are often performed in modifier-free conditions.

1. Introduction

Supercritical solvents and, in particular, supercritical carbon dioxide (sc-CO₂) has gained world-wide attention in recent decades as a green processing media for the impregnation of porous matrices like wood, textile, monolithic and porous polymers with various solutes such as pharmaceutically active compounds, pesticides, or organometallic compounds, among others [1,2,3]. Sc-CO₂ is an environmentally friendly alternative to conventional organic solvents. The reasons for which sc-CO₂ has attracted attention are as follows: (1) its high diffusivity and low viscosity diminish mass transfer limitations; (2) mild conditions (above 31 °C, 73 bar), that are easily achievable with standard lab and industrial compressors and fittings, allow the processing of sensitive and thermolabile materials; (3) non-flammability, low price, and GRAS status (generally recognized as safe) by the European Food Safety Authority and Food and Drug Administration (USA). Furthermore, impregnated materials can be recovered free of solvent residues by simple depressurization. As a result of all these advantages, sc-CO₂-based processes have a great potential, particularly, for impregnation in the biomedical field [4,5,6,7,8]. Two mechanisms play the major roles in the impregnation processes: adsorption of the solutes on the matrix surface and precipitation upon solvent evaporation [1]. In this work, we mainly focus on the first mechanism and aim to quantify adsorption interactions between a given solute on a given matrix, which contributes to the rational development of such impregnation processes.

Broadly, methods for studying adsorption equilibria in supercritical solvents can be classified into static and dynamic ones [9]. Although static methods offer superior precision, they are often too time consuming for screening of the affinity of numerous solutes towards a selected matrix. To ensure a high throughput, dynamic methods are employed instead [10,11,12,13].

Dynamic methods rely on the fundamental link between thermodynamic and kinetic properties. One common example is the partitioning of the solute between an adsorbent and flowing fluid in chromatographic experiments, where retention factor $k_0$ can be used to obtain distribution coefficients and adsorption isotherms [14].

Thus, supercritical fluid chromatography (SFC) using a porous matrix of interest as stationary phase and solutes injected into pure sc-CO₂ (mobile phase) can provide direct access to adsorption equilibria. Such an implementation of SFC is referred to as inverse SFC (iSFC), stressing the fact that retention characteristics of a single solute are studied. Using fundamental thermodynamic relations for the retention factor $k_0$ (Subscript “0” refers to the condition of using pure sc-CO₂ as mobile phase) [14,15], iSFC with pure sc-CO₂ as a mobile phase has been used to determine distribution coefficients, adsorption enthalpies, and partial molar volumes at infinite dilution on activated carbons, resins, and plain and modified silicas (see [14]).

Typically, even moderately polar solutes demonstrate a high affinity to porous matrices [10], so their retention time is extremely large and the peak broadening is severe with pure sc-CO₂, given its comparably low solvation power. Thus, it is not possible to determine the precise retention factor. Addition of a few percent of organic solvents, so-called modifiers, are commonly employed to remedy this issue. This—at first sight, just a practical difficulty—has serious negative consequences for the interpretation of the retention factor in terms of thermodynamic quantities: the retention factor $k$ measured in systems with modifier cannot be directly related to the retention factor $k_0$ in modifier-free conditions.

Our aim is to develop models for the estimation of retention factors for polar solutes in pure sc-CO₂ using retention data obtained in sc-CO₂/modifier mixtures, that is, to deduce $k_0$ values from experimentally measurable retention factors, $k$. The developed models were validated against experimental data.

2. Materials and Methods

2.1. Stationary Phase, Mobile Phase, and Solutes

Four silica-based porous materials were employed as stationary phases in this work (

Table 1. Textural properties of the porous materials employed as stationary phases.

Column	Kromasil® 60-5-SIL	Kromasil® 100-5-SIL	Kromasil® 300-5-SIL	Aerogel
Abbreviation	SIL-60	SIL-100	SIL-300	SIL-Aerogel
* Particle size, μm	6.2	6.1	6.2	17.4
** Mean pore diameter, nm	4	6	23	27
** Specific surface area, m2 g−1	540	320	110	858
** Specific pore volume, cm3 g−1	1.1	0.9	0.9	5.2

* Measured by electrical sensing zone particle analyzer and scanning electron microscope. ** Measured by N₂-adsorption/desorption with Nova 3000e surface area and pore analyzer (Quantachrome instruments); BET and BJH methods were employed for surface area and pore volume, respectively.

). Three commercial silica gels, Kromasil^® 60-5-SIL, 100-5-SIL, and 300-5-SIL, were supplied by AkzoNobel (Bohus, Sweden). They are referred to as “SIL-60”, “SIL-100”, and “SIL-300”, respectively, in the following. The fourth material is a silica aerogel produced in our laboratory, referred to as “SIL-Aerogel”. The porous materials were packed in commercial HPLC columns as reported in our previous work [10]. A mixture of supercritical CO₂ (Westfalen AG, Münster, Germany) and methanol (Roth, Karlsruhe, Germany) was employed as mobile phase. Phenol (Fluka, Taufkirchen, Germany), toluene (Carl Roth, Karlsruhe, Germany), naproxen, nicotinamide, ketoprofen (Merck, Kenilworth, NJ, USA), and rosmarinic acid (Sigma-Aldric, Steinheim, Germany) were employed as solutes. All solutes were dissolved in methanol to a concentration of 0.1 g L⁻¹. The hold-up time $t_M$ was measured with a solution of nitrous oxide (N₂O) (Westfalen AG, Münster, Germany) in methanol. The solution was prepared by bubbling N₂O through a metal solvent filter in methanol for at least five minutes.

2.2. Chromatographic System and Measurements

SFC system: both retention and adsorption isotherm measurements were carried out with Waters Acquity UPC² supercritical fluid chromatography system equipped with a photodiode array (PDA) detector (Waters Corporation, Milford, CT, USA).

Hold-up time and retention data: Mobile phase was pumped at a flow rate of 2 mL min⁻¹ with 0 to 30 vol.% methanol. Columns packed with porous materials were fixed in an insulated container equipped with a preheater enabling temperature control in increments of 0.1 ° C. The pressure was controlled by back pressure regulators with an accuracy of 0.5 bar. The samples were injected into the SFC system and had volumes of 2 µL. The PDA detector signal was recorded for each compound at its specific working wavelength: phenol at 210 nm, toluene 240 nm, naproxen at 220 nm, nicotinamide at 210 nm, ketoprofen at 254 nm, rosmarinic acid at 330 nm, and nitrous oxide at 195 nm [16]. Empower^® 3 software [17] was used to operate the SFC system.

The time at which the maximum signal appears was taken as the retention time, $t_R$. The retention time of nitrous oxide was employed as the hold-up time, $t_M$. Retention and hold-up times of all solutes were measured at least in triplicate and averaged. The averaged values were used to calculate the retention factor, as in Equation (1):

$$k=\frac{t_R-t_M}{t_M}$$

(1)

The methanol concentration was varied from 0 to 30 vol.%. Retention factors were plotted against methanol concentration to generate a retention profile.

Determination of methanol adsorption isotherms: Adsorption isotherms of methanol on all porous matrices in sc-CO₂ were measured by frontal analysis chromatography with the staircase method. This method was used in several previous studies [11,18].

Breakthrough curves were measured in the SFC systems for a series of methanol concentrations for each silica material. By integrating the area above the breakthrough curves, the maximum loadings at specific methanol concentrations were determined. Obtained loadings were fitted to the Langmuir equation to estimate the equilibrium constant $K_{eq}$. The complete experimental procedure and integration method are detailed in Supplementary Materials (Section S1, Figure S1, Table S1).

3. Retention Modeling

3.1. Model Development

For a retention model to be useful, it should relate the retention factor $k$ to operating conditions such as pressure, temperature, and composition of the mobile phase. Once such a model is obtained, the retention factor can then be related to the distribution constant of a solute between mobile phase and stationary phases (Equation (2)) and other thermodynamic quantities [14].

$$k=D\beta =\frac{c_s}{c_m}\beta$$

(2)

Here, $D$ is the distribution coefficient, $c_s$ and $c_m$ are the concentrations of a solute (mg mL⁻¹) in the stationary and mobile phases, and $\beta$ is the ratio of stationary to mobile phase volumes [19].

With a retention model in hand which is valid at zero modifier concentration, the retention factor $k_0$ can be then extracted from experimental retention data. In this work, two classes of retention models were developed based on different assumptions.

3.1.1. Dual Influential Model

The first category of models is based on the idea that modifier influences the solute retention both in the mobile and stationary phases. We assume that the ratio of the solute concentration in the mobile phase with modifier $c_m$ over the corresponding concentration without modifier is proportional to the modifier concentration $c_{mod}$ (Equation (3)). Likewise, the solute concentration in the stationary phase is increased proportional to the modifier loading in the stationary phase (Equation (4)). The loading is assumed to be independent of the presence of the solute and is a function of the pressure and temperature of the mobile phase. This class of models is referred to as dual influential models (DIM). Equations (3) and (4) are suggested in the spirit of the Setschenow equation for the salting out effect and the log-linear model for solubility in mixed solvents [20]:

$$\lg \frac{c_m}{c_{m,0}}=a{c}_{mod}$$

(3)

$$\lg \frac{c_s}{c_{s,0}}=b{q}_{mod}$$

(4)

where $c_{m,0}$ and $c_{s,0}$ are the concentrations of a solute in mobile and stationary phases at zero modifier concentration, $c_{mod}$ is the modifier concentrations (mg mL⁻¹) in mobile phase, $q_{mod}$ is the loading (mg mL⁻¹) of modifier on stationary phase by adsorption, and both $a$ and $b$ are empirical parameters (mL g⁻¹) that can be regarded as a strength of the modifier effect on the mobile and stationary phases, respectively. When the absolute values of $a$ and $b$ are large, the solute concentrations are more sensitive to small changes of modifier concentration. The values of $a$ and $b$ depend on the nature of the solute, the mobile and stationary phase, and the operational conditions (temperature, pressure).

If we now assume the simplest case, namely a linear isotherm for the modifier adsorption equilibrium, we get Equation (5), as follows:

$$q_{mod}=K_{eq, lin}{c}_{mod}$$

(5)

where $K_{eq, lin}$ is the equilibrium constant. Equations (2)–(5) can be rearranged into the following expression for the retention factors $k$ at the modifier concentration $c_{mod}$ and $k_0$ is the retention factor at $c_{mod}=0$:

$$\lg k=\lg k_0-A{c}_{mod}$$

(6)

where $A=a-b{K}_{eq,lin}$. The constant $A$ is the overall strength of the modifier effect, which accounts for the change in the retention factor from its value $k_0$ under modifier-free conditions.

Another natural extension is to plug a non-linear adsorption isotherm into the retention model. Several previous works have demonstrated that the adsorption of individual organic solvents in supercritical CO₂ follows Langmuir-type isotherm models [13,21,22,23,24].

In this study, we also demonstrate that Langmuir isotherm, Equation (7), provides an acceptable fit to the experimental data for methanol on silica-based materials (see Section S1 in Supplementary Materials).

$$q_{mod}=\frac{q_{max}{K}_{eq,L}{c}_{mod}}{1+K_{eq,L}{c}_{mod}}$$

(7)

The fractional surface coverage $\theta$ of the stationary phase by modifier molecules is given as follows:

$$\theta =\frac{q_{mod}}{q_{max}}=\frac{K_{eq,L}{c}_{mod}}{1+K_{eq,L}{c}_{mod}}$$

(8)

Here, $K_{eq,L}$ is the Langmuir equilibrium constant (mg mL⁻¹) and $q_{max}$ is the maximum modifier loading (mg mL⁻¹) at monolayer coverage. Several $K_{eq,L}$ values were measured as a part of our preliminary work (see Section S1 in Supplementary Materials).

Equation (8) together with Equations (2)–(4) results in the retention model given by Equation (9).

$$\lg k=\lg k_0+B\theta -a{c}_{mod}$$

(9)

$$\lg k=\lg k_0+B\frac{K_{eq, L}{c}_{mod}}{1+K_{eq, L}{c}_{mod}}-a{c}_{mod}$$

(10)

where $B=b{q}_{max}$. Equation (10) gives a concentration-explicit form of Equation (9) and will be tested for the ability to estimate $k_0$ values along with linear-DIM model, Equation (6).

3.1.2. Mixed-Retention Model

Next, we develop and study another category of retention models which are based on the following assumption: the solute retention is due to independent adsorption on two types of active sites of the stationary phase, free adsorption sites as well as the active sites covered by adsorbed modifier (Figure 1). We call such models mixed retention models (MRM). Experimentally observed retention is a joint effect of the retentions due to both sites.

For plain silica, the adsorption sites are silanol groups and methanol molecules adsorbed on the silanols. Since the fraction of the methanol-covered sites is identical to the modifier surface coverage, $\theta$, the retention factor can be expressed by a weighted sum of the retention factor at zero modifier concentration, $k_0$, and the retention factor for completely covered surface, $k_c$ (subscript $c$ stands for covered):

$$k=k_0\left(1-\theta \right)+k_c\theta$$

(11)

If we rearrange Equation (11), a linear relation between $k$ and $\theta$ is seen:

$$k=\left(k_c-k_0\right)\theta +k_0$$

(12)

However, according to the collected experimental data of $k$ and $\theta$ in this work, the linear relationship does not hold for the whole range of modifier concentrations. Details will be given later in Section 4.2. Testing various modifications of Equation (12) that better fit the experimental data, we introduced an empirical parameter $M$ to outweigh the influence of covered active sites; by introducing the term $M{c}_{mod}$ to the denominator, the influence of modifier on retention factor reduces at higher modifier concentration. The extended MRM is given by the following:

$$k=\frac{k_0\left(1-\theta \right)+k_c\left(1+M\right)\theta }{1+M{c}_{mod}}$$

(13)

The terms $\left(1+M\right)$ in the numerator and $M$ in the denominator are introduced to make the original model by Equation (11) a special case of the extended model if $M=0$.

The MRM given by Equation (13) can also be written in a concentration-explicit form by substitution of a certain adsorption isotherm $\theta =f\left(c_{mod}\right)$. For example, for the Langmuir isotherm, we obtain:

$$k=\frac{k_0+K_{eq,L}{k}_c\left(1+M\right)c_{mod}}{\left(1+K_{eq,L}{c}_{mod}\right)\left(1+M{c}_{mod}\right)}$$

(14)

3.2. Evaluation of the Models

As explained above, we have discussed two categories of retention models, DIMs and MRMs. Each category contains several models, derived with different isotherm types or extended versions. Furthermore, even with the same model, fits can be performed using different input data arrangements and lead to significantly different results. To keep the text clearer, we introduce the following notation. The name of each method contains three parts: (1) either a DIM or MRM, (2) either linear (Lin) or Langmuir isotherm (Lan), and (3) the symbols of input data. Taken together, all these elements form what we call a method.

Table 2. Description of employed methods: each method is a combination of a model (DIM or MRM), adsorption isotherm, and a set of fitted parameters.

Method	Equation Used	Input Data	Fitted Parameters
DIM-Lin-kc	$\lg k=\lg k_0-A{c}_{mod}$	$k,c_{mod}$	$A$,$k_0$
DIM-Lan-kcK	$\lg k=\lg k_0+B\frac{K_{eq,L}{c}_{mod}}{1+K_{eq,L}{c}_{mod}}-a{c}_{mod}$	$k, c_{mod}, K_{eq,L}$	$B$, $a$, $k_0$
DIM-Lan-kc	$\lg k=\lg k_0+B\frac{K_{eq,L}{c}_{mod}}{1+K_{eq,L}{c}_{mod}}-a{c}_{mod}$	$k,c_{mod}$	$B$, $a$, $k_0$, $K_{eq,L}$
DIM-Lan-kcK’	$\lg k=\lg k_0+B\frac{K_{eq,L}{c}_{mod}}{1+K_{eq,L}{c}_{mod}}-a{c}_{mod}$	$k, c_{mod}, K_{eq, L}^{corr}$	$B$, $a$, $k_0$
MRM-Lan-kcK	$k=\frac{k_0+K_{eq,L}{k}_c\left(1+M\right)c_{mod}}{\left(1+K_{eq,L}{c}_{mod}\right)\left(1+M{c}_{mod}\right)}$	$k, c_{mod}, K_{eq,L}$	$k_0$, $k_c$, $M$
MRM-Lan-kcKM	$k=\frac{k_0+K_{eq,L}{k}_c\left(1+M\right)c_{mod}}{\left(1+K_{eq,L}\bullet c_{mod}\right)\left(1+M{c}_{mod}\right)}$	$k, c_{mod}, K_{eq,L},M$	$k_0$, $k_c$

gives an overview of the employed methods. All regressions (linear and nonlinear) were performed in Origin 8.5.1 [25], through chi-square minimization with a tolerance of 10⁻⁹, and a statistical weighting method was applied. A summary of fitting parameters applied for each method and the ranges of R² values of regressions are given in

Table 3. Comparison of methods based on DIMs an MRMs.

Method	Number of Fitting Parameters	Range of R2	Experimental Efforts
DIM-Lin-kc	2 ($A$ and $k_0$)	0.74–0.94	Low
DIM-Lan-kcK	3 ($B$, $\alpha$, and $k_0$)	0.83–1.00	High
DIM-Lan-kc	4 ($B$, $\alpha$, $k_0$ and $K_{eq,L}$)	0.87–1.00	Low
DIM-Lan-kcK’	3 ($B$, $\alpha$, and $k_0$)	0.84–1.00	Low
MRM-Lan-kcK	3 ($k_0$, $k_c$, $M$)	(0.86) *–1.00	High
MRM-Lan-kcKM	2 ($k_0$, $k_c$)	0.69–1.00	Low

* 0.86 is the lowest R² among those with successful fitting. However, several fittings failed by this method.

.

4. Results and Discussion

4.1. Applying Dual Influential Models

In this section, the DIMs based on linear and Langmuir isotherms were compared by analyzing the results of the fittings to experimental data. The appropriate method to proceed with the regression analysis according to different situations will be suggested after the comparison and discussion.

4.1.1. Model Assuming Linear Isotherms for Modifier Adsorption

Method DIM-Lin-kc applied Equation (6) for data fitting. Equation (6) was derived based on the assumption that the loading of modifier in the stationary phase is proportional to the modifier concentration in the mobile phase. This assumption was applied as a first attempt, intended to obtain the simplest model with the least number of fitting parameters, although it is known that methanol adsorption on porous silica in sc-CO₂ does not follow a linear isotherm across the whole range of methanol fraction [26].

The results of the fitting for retention factors of phenol on silicas are demonstrated in Figure 2, and those for other solutes on SIL-60 are in Figure 3. The model fittings were all successful in the sense that either the chi-square minimization reached the tolerance of 10⁻⁹ or the chi-square no longer changed. The adjusted R² (adj. R²) for all cases were in the range of 0.74–0.94, and the relative standard errors (RSEs) of the fitting parameter $A$ and $k_0$ were in the range of 7–16%. However, it can be observed from the figures that the model fittings diverge from the experimental data at low concentrations of methanol. This leads to the underestimations of the $k_0$ values (y-intercept) in all cases regardless of physical properties of stationary phase (Figure 2) and the solute nature (Figure 3). These observations point to inadequacy of the model assumptions: as expected, the use of linear isotherms to describe the modifier adsorption on porous silica materials is an oversimplification and cannot be recommended for estimation of $k_0$.

Note that Equation (6), derived in this work, has a certain similarity to the common retention model for reversed phase chromatography, where the logarithm of the retention factor of a solute is linearly related to the organic solvent fraction and the intercept is the logarithm of the retention factor in pure water as mobile phase [27]. This simple model was found to be not accurate enough for certain systems. Efforts were made to improve the model by introducing quadratic terms to the right-hand side [28]. Another natural way, however, is to apply non-linear adsorption isotherms.

4.1.2. Models Assuming Langmuir Isotherms for Modifier Adsorption

As mentioned in Section 3.1.1, Langmuir isotherm was adopted for constructing retention models to describe the modifier adsorption in this work. Three methods based on Langmuir isotherm to analyze the retention data were established. Together with the linear-isotherm-based method discussed in the previous section, in total four methods based on DIMs were studied. Estimated $k_0$ and adj. R² for the four methods are listed in

Table 4. Regression results for DIMs.

Stationary Phase	Solute	DIM-Lin-kc		DIM-Lan-kcK			DIM-Lan-kc			DIM-Lan-kcK’
		${\mathit{k}}_0$	R2	${\mathit{k}}_0$	${\mathit{K}}_{\mathit{e}\mathit{q},\mathit{L}}$	R2	${\mathit{k}}_0$	${\mathit{K}}_{\mathit{e}\mathit{q},\mathit{L}}$	R2	${\mathit{k}}_0$	${\mathit{K}}_{\mathit{e}\mathit{q}, \mathit{L}}^{\mathit{c}\mathit{o}\mathit{r}\mathit{r}}$	R2	${\mathit{k}}_0^{\left(\mathit{D}\mathit{I}\mathit{M}-\mathit{L}\mathit{a}\mathit{n}-\mathit{k}\mathit{c}\mathit{K}’\right)}/{\mathit{k}}_0^{\left(\mathit{D}\mathit{I}\mathit{M}-\mathit{L}\mathit{a}\mathit{n}-\mathit{k}\mathit{c}\mathit{K}\right)}$
SIL-60	Phenol	4.4 ± 0.3	0.92	6.2 ± 0.2	16.5	0.99	8.1 ± 0.1	62 ± 2	1.00	5.9 ± 0.2	11.7	0.99	0.96
	Naproxen	25 ± 5	0.74	36 ± 4		0.93	330 ± 19	170 ± 4	1.00	32 ± 4		0.91	0.91
	Nicotinamide	42 ± 5	0.92	304 ± 16		1.00	(1.4 ± 0.2) × 103	40 ± 3	1.00	219 ± 14		1.00	0.72
	Rosmarinic acid	483 ± 81	0.94	(1.02 ± 0.08) × 104		1.00	(2.9 ± 1.2) × 105	45 ± 4	1.00	(5.6 ± 0.4) × 103		1.00	0.55
	Toluene	0.2	0.89	0.2		0.97	0.2	308 ± 131	0.98	0.2		0.97	1.00
SIL-100	Phenol	3.0 ± 0.2	0.89	4.1 ± 0.1	11.2	0.99	5.6 ± 0.1	58 ± 2	1.00	4.1 ± 0.1	11.1	0.99	1.01
SIL-300	Phenol	0.9 ± 0.1	0.78	1.1	5.3	0.97	1.4	24 ± 3	1.00	1.1	4.6	0.97	0.99
SIL-Aerogel	Phenol	2.0 ± 0.2	0.90	2.6 ± 0.1	7.1	0.98	3.9 ± 0.1	63 ± 3	1.00	2.3 ± 0.1	11.9	0.96	0.90
	Naproxen	96 ± 23	0.93	139 ± 8		1.00	641 ± 178	26 ± 3	1.00	185 ± 10		1.00	1.33
	Toluene	0.1	0.74	0.1		0.83	0.1	433 ± 487	0.87	0.1		0.84	1.05

. All experimental data were obtained at 40 °C, 200 bar. The method DIM-Lan-kcK’ employs a corrected $K_{eq, L}^{corr}$ as input instead of measured $K_{eq, L}$. This method was established after analyzing the results from the first three methods and is described in the following.

Regarding the aim of selecting the best method, three aspects were considered: (1) experimental efforts, (2) goodness of fitting, and (3) number of fitting parameters or problem of overfitting.

The DIM-Lan-kcK method relies on experimentally measured $K_{eq, L}$ as input data for fitting; thus, we kept a small number of fitting parameters to avoid overfitting. However, measurements of $K_{eq, L}$ require intensive efforts considering both time and equipment upgrades [13,22]. In the method DIM-Lan-kc, we attempted to avoid such troublesome experimental measurements by making $K_{eq, L}$ a fitting parameter. The fitted $K_{eq, L}$ values were found to be significantly larger than the measured ones. In addition, the equilibrium constants should remain nearly constant for each stationary phase at a constant temperature regardless of the used solutes. Contrary to this, the fitted values of $K_{eq, L}$ obtained for different solutes deviate largely for SIL-60 and SIL-aerogel (data not shown).

We surmise that overfitting of the retention data is a possible reason. Nevertheless, we observed that the fitted $K_{eq, L}$ values obtained in the DIM-Lan-kc method for phenol correlate well with the experimentally measured $K_{eq, L}$ values. Thus, the $K_{eq, L}$ obtained from the DIM-Lan-kc method should be first corrected by multiplying by an empirical factor, $F$. We call this value $K_{eq, L}^{corr}$:

$$K_{eq, L}^{corr}=F{K}_{eq, L}^{’}$$

(15)

In the DIM-Lan-kcK’ method, the corrected adsorption constant from Equation (15) was employed to solve the overfitting problem of the DIM-Lan-kc method. An estimation for the correction factor $F$ was made through a paired t-test (see details in Supplementary Materials, Section S2, Tables S2 and S3) and found to be 0.19. The value of $k_0$ estimated by the DIM-Lan-kcK’ method is close (55–133%, Table 4) to that from the method DIM-Lan-kcK, where experimental values of $K_{eq, L}$ are used.

Compared to linear adsorption isotherm, the use of Langmuir isotherms showed a significantly better correlation, as expected, due to a greater number of fitting parameters in the latter case. However, although the DIM-Lan-kc method had one more fitting parameter than DIM-Lan-kcK and DIM-Lan-kcK’, the adj. R² did not show a visible increase. Another observation regarding the goodness of fitting is that the regression for toluene always results in the smallest adj. R². Toluene was selected as it is the only solute for which $k_0$ could be measured experimentally due to its low affinity to silica. However, because of the same reason, retention times of toluene are close to the hold-up time, resulting in large uncertainties of the retention factor. Therefore, the low adj. R² for toluene should not be interpreted as a failure of the model.

To summarize, the DIM-Lin-kc method is not recommended in any case because it has the worst goodness of fitting. Having the same level of goodness of fitting, the DIM-Lan-kcK’ method shows significant advantages compared to DIM-Lan-kcK for less experimental efforts and to DIM-Lan-kc for less risk of overfitting (Figure 4 and Figure 5).

Based on the above results, several issues shall be addressed when applying the DIMs for estimating $k_0$ from retention data in the presence of a modifier. First of all, when the methanol adsorption isotherm on the to-be-tested stationary phase in supercritical CO₂ is available at the required temperature and pressure, the method DIM-Lan-kcK can be applied with the measured $K_{eq,L}$ for the most reliable estimation. However, in most cases, adsorption isotherms are not available. If so, phenol should be injected into the column as a model solute at various methanol concentration. The $K_{eq,L}$ can be obtained afterwards by regression using DIM-Lan-kc on the retention data of phenol. The value of A $K_{eq, L}^{corr}$ can be then calculated by multiplying by the correction factor $F$. DIM-Lan-kcK’ can be finally applied for estimation of $k_0$ of any solutes on the stationary phase of this column using the $K_{eq, L}^{corr}$.

4.2. Applying Mixed-Retention Models

According to the original MRM, Equation (12), the retention factor $k$ varies linearly with surface coverage $\theta$:

$$k=\left(k_c-k_0\right)\theta +k_0$$

(16)

As seen in the $k$ vs. $\theta$ plot (Figure 6), the experimental data does not show linear correlations among any solute/stationary phase pairs. Only the data points of each solute/stationary phase pair at relatively lower modifier concentrations can be treated as a linear function. The reason might be that the assumption of single layer adsorption to apply Langmuir isotherm does not hold true. Studies have demonstrated the adsorption of alcohols on silica materials could exhibit multilayer behavior above a certain concentration [26,29,30]. If this is the case, two issues should be considered. Firstly, the calculated $\theta$ values are only hypothetical and do not have any physical meaning. Secondly, the proposed model, Equation (16), is based on the theory that when methanol is adsorbed as a single layer, the solute molecules are attracted to the methanol molecules, which leads to $k_c$. However, when multilayer adsorption appears, the affinity between the solute and adsorbed methanol changes, which leads to a different $k_c^{’}$. In this case, $k_c$ depends on the modifier concentration, and the model assumptions break down. Because of these reasons, an advanced model Equation (14) with an empirical factor $M$ was suggested. As discussed above for the DIMs, when the experimental value of $K_{eq,L}$ is not available, an estimated value can be taken. The results from regression by MRM-Lan-kcK and MRM-Lan-kcKM methods are given in

Table 5. Regression results for MRMs.

Stationary Phase	Solute	MRM-Lan-kcK				MRM-Lan-kcKM
		${\mathit{k}}_0$	$\mathit{M}$	${\mathit{K}}_{\mathit{e}\mathit{q},\mathit{L}}$	R2	${\mathit{k}}_0$	$\mathit{M}$	${\mathit{K}}_{\mathit{e}\mathit{q},\mathit{L}}$	R2	${\mathit{k}}_0^{\left(\mathit{M}\mathit{R}\mathit{M}-\mathit{L}\mathit{a}\mathit{n}-\mathit{k}\mathit{c}\mathit{K}\mathit{M}\right)}/{\mathit{k}}_0^{\left(\mathit{D}\mathit{I}\mathit{M}-\mathit{L}\mathit{a}\mathit{n}-\mathit{k}\mathit{c}\mathit{K}\right)}$
SIL-60	Phenol	7.6	113 ± 13	16.5	1.00	8.3 ± 0.1	150.0	20.0	1.00	1.35
	Naproxen	-	-		-	31.7			0.90	0.89
	Nicotinamide	-	-		-	277 ± 10			0.99	0.91
	Rosmarinic acid	-	-		-	915 ± 276			0.69	0.09
	Toluene	0.2	184 ± 41		0.98	0.2			0.98	1.03
SIL-100	Phenol	5.4 ± 0.2	135 ± 9	11.2	1.00	5.6 ± 0.1	150.0	20.0	1.00	1.39
SIL-300	Phenol	1.7 ± 0.2	191 ± 42	5.3	0.98	1.5 ± 0.1	150.0	20.0	0.97	1.35
SIL-Aerogel	Phenol	3.8 ± 0.1	152 ± 10	7.1	1.00	3.78 ± 0.04	150.0	20.0	1.00	1.47
	Naproxen	-	-		-	155 ± 11			0.97	1.12
	Toluene	0.09 ± 0.01	200 ± 139		0.86	0.09 ± 0.01			0.86	0.99

.

We have observed that when Equation (14) was applied using the experimental measured $K_{eq,L}$ (MRM-Lan-kcK), the data could not be regressed and the parameters did not converge. The suspected reason was over-parameterization, which leads to deterioration of estimation accuracy, and the resulted estimations and errors were unreliable. However, the fitting worked in all cases for phenol and toluene as solutes and the estimated $k_0$ values were always comparable to those from DIM-Lan-kcK. The parameter $M$ was found to lie in the range 100–200.

Although the retention factors of toluene and phenol are very different at zero modifier concentration, they have similar $M$ values. We may surmise that $M$ is not solute-dependent, but modifier-dependent. This is reasonable because $M$ was introduced to correct the strength of the modifier influence at high modifier concentration where multilayer adsorption may take place and the influence of solutes shall be minimal. Therefore, to solve the overparameterization problem, $M$ was fixed to 150 in all cases for the next method, MRM-Lan-kcKM.

Unlike developing method DIM-Lan-kcK’ based on the result from DIM-Lan-kc, it is not possible to treat the $K_{eq,L}$ as a fitting parameter in the MRMs: it is not reasonable to add an extra fitting parameter to a model which is already overparameterized. We suggest a rough estimation of $K_{eq,L}$ = 20 mL g⁻¹ as a representative value across all studied silica materials. With these estimates for $M$ and $K_{eq,L}$, the regressions were performed and the results show reasonable agreement (89–147%) with the method DIM-Lan-kcK (Table 5). One exception was the case of rosmarinic acid, which has an exceptionally high affinity to silica. In this case, the estimations of $k_0$ by MRM are not reliable. For the second strongest retained solute, nicotinamide, the DIM-Lan and MRM-Lan methods give consistent estimations (304 ± 16 and 277 ± 10, respectively), even when arbitrary selected parameters $K_{eq,L}$ and $M$ were taken.

To summarize, if $K_{eq,L}$ is not available, the method MRM-Lan-kcKM is recommended; it provides satisfying estimations close to DIM-Lan-kcK without over-parameterization.

4.3. Case Study: Temperature Influence on Retention

To validate the methods explained in the previous chapters and to study the temperature influence on the retention factors, injections of ketoprofen were performed at constant CO₂ density (807 g L⁻¹). Four temperatures were investigated and the back pressure of the system was adjusted accordingly. The concentration of the modifier was very low and considered as having no influence on fluid density. The experimental data and the model fitting with the MRM-Lan-kcKM method is given in Figure 7, while the estimated $k_0$ values using the methods DIM-Lan-kcK, DIM-Lan-kcK’ and MRM-Lan-kcKM are presented in Figure 8.

The temperature dependence of adsorption in supercritical CO₂ has been studied in many previous works [31,32,33,34,35]. However, all of these works were done by varying temperatures at a fixed pressure. No matter which solutes or stationary phases were studied, the same conclusion was drawn: that the adsorbed amount increases with increasing temperature. This phenomenon was explained mainly by the density influence. With increasing temperature, the density and solvation power of CO₂ is lowered, favoring solute adsorption on the solid surface. In the current work, we eliminated this cause by fixing the density to a constant value, and we found that the $k_0$ decreases with increasing temperature. When temperature increases, the vapor pressure of ketoprofen increases, and therefore, ketoprofen shows lower affinity to the stationary phase at higher temperatures. This behavior follows the general temperature dependence of adsorption in gases and liquids. It is also seen in Figure 7 that the four lines come closer at higher $c_{mod}$, pointing to the fact that temperature has less influence at higher $c_{mod}$ on the adsorption of solute to SIL-Aerogel. One explanation is that the higher modifier content leads to higher solubility of solutes in the mobile phase, which dominates over the temperature effect. Another explanation could be that the temperature influence on the adsorption of modifier and solute counteract each other: at lower temperatures, the adsorption of both modifier and solute is higher. At the same time, since more modifier is adsorbed, it is covering a larger surface and inhibiting the adsorption of the solute. This counteraction is more obvious at higher methanol concentrations.

Observed in Figure 8, as expected, the estimated $k_0$ values from all three methods were constant taking into account the margin of error. This result validated the two methods DIM-Lan-kcK’ and MRM-Lan-kcKM, which do not require prior determinations of adsorption isotherms of modifiers. A quick summary of the key estimations of the two methods: Both methods assume the modifier adsorption follows Langmuir isotherm. To apply DIM-Lan-kcK’, the $K_{eq, L}^{corr}$ can be calculated from the retention data of phenol on a selected column. A correction factor $F=0.19$ shall be adopted. The obtained $K_{eq, L}^{corr}$ can be then used further for other solutes. To apply MRM-Lan-kcKM, rough estimations of $K_{eq,L}$ = 20 mL g⁻¹ and $M=150$ should be taken.

5. Conclusions

This work reports on experimental retention factors for moderately polar solutes on four silica-based porous matrices used as stationary phases in packed-column SFC. Low polarity of pure sc-CO₂ makes elution of the solutes only possible if modifier (methanol) is added to the mobile phase. To extrapolate the retention factors into experimentally inaccessible modifier-free conditions, two classes of empirical retention models were developed. The first class relies on log-linear relations for the equilibrium concentration of the solute in both mobile and stationary phases. The second class interprets the overall retention as a sum result of the retention on native and covered-by-modifier adsorption sites. To obtain analytical expressions for both classes, a functional form for adsorption isotherms of modifier on the porous matrix is required. We demonstrated that linear adsorption isotherm is insufficient to describe experimental results. Langmuir isotherm yields qualitatively correct behavior. To obtain quantitative agreement and to avoid overfitting of the models, we suggested how to estimate the adsorption coefficient of the modifier. An additional advantage is that time-consuming prior determinations of adsorption isotherms of modifier can be avoided. We apply the developed models to the adsorption of a model drug, ketoprofen, on silica aerogel and deduce the retention factor in modifier-free conditions based on retention data obtained with modifier. Since the impregnation processes are commonly performed in modifier-free conditions, the developed models may facilitate the screening of porous matrices with a required affinity towards solutes of interest.