Predictions of overloaded concentration proﬁles in supercritical ﬂuid chromatography

Here, overloaded concentration proﬁles were predicted in supercritical ﬂuid chromatography using a combined two-dimensional heat and mass transfer model. The heat balance equation provided the temperature and pressure proﬁles inside the column. From this the density, viscosity, and mobile phase velocity proﬁles in the column were calculated. The adsorption model is here expressed as a function of the density and temperature of the mobile phase. The model system consisted of a Kromasil Diol column packed with 2.2-μm particles (i.e., a UHPSFC column) and the solute was phenol eluted with neat carbon dioxide at three different outlet pressures and ﬁve different mobile phase ﬂow rates. The proposed model successfully predicted the eluted concentration proﬁles in all experimental runs with good agreement even with high-density drops along the column. It could be concluded that the radial temperature and density gradients did not signiﬁcantly inﬂuence the overloaded concentration elution proﬁles.


Introduction
Interest in supercritical fluid chromatography (SFC) is still growing, mainly due to the method's potential for environmental friendliness as well as fast and efficient separations.The diffusion coefficients in SFC can in some cases be several times greater than those observed in HPLC.Generally, the viscosity of the mobile phase in SFC is much lower than that in LC and, as a result, we observe a much lower pressure drop along the column [1] .
However, SFC is more complex than liquid chromatography due to the high compressibility of the supercritical fluid [ 1 , 2 ].The isenthalpic expansion of the mobile phase leads to column cooling, and significant temperature gradients can form in the axial and radial directions.These gradients influence the density, viscosity, and velocity of the mobile phase as well as the diffusion and dispersion coefficients, which cannot be assumed to be constant as is the case in HPLC.With the very high expansion of the mobile phase, phenomena similar to those in ultra-high-pressure liquid chromatography (UHPLC) can be observed, but their complexity is much greater.It is therefore difficult to fully understand all the loaded concentration profiles by assigning them to the section of the column with the same pressure as the average pressure applied during their determination.A linear change in pressure along the column was assumed.The velocity of the mobile phase in each section of the column was calculated, additionally assuming a linear change in the temperature of the mobile phase.
Wenda and Rajendran [7] considered the separation of the enantiomers of flurbiprofen on the chiral stationary phase in SFC under both linear and nonlinear conditions.They determined the competitive Langmuir and Bi-Langmuir isotherm models under different experimental conditions using the inverse method (only the mass transfer model was used).The experiments were conducted in such a way that the pressure drop did not exceed 35 bar while the maximum density drop was 15 g L -1 .A procedure similar to that of Wenda and Rajendran has been used by others [8][9][10][11][12] .In this procedure, the overloaded bands are modeled but the variations in the isotherm parameters and the mobile phase velocity with temperature and pressure are neglected or only the velocity changes in the axial direction are included.As far as we know, no prediction has been made of the overloaded profiles in SFC with large temperature and density changes included in the calculation of the heat transport model.
The overarching aim of this study is to predict overloaded concentration profiles under non-isopycnic conditions in SFC and, when doing so, to account for all changes in the mobile phase properties in the axial and radial directions.The overloaded concentration profiles will be calculated using the two-dimensional equilibrium dispersive model.Here an adsorption model, which is a function of the temperature and density of the mobile phase, must be applied.The mathematical model, as well as the means of estimating the parameters and calculating the temperature distribution and overloaded concentration profiles, will be validated using basic experimental data obtained at different mobile phase flow rates and outlet pressures, with the significant changes in temperature and especially density being observed.A further aim is to examine the thermal effects occurring in the SFC column and their possible impact on the studied model cases.

Mathematical model
Proper modeling of the concentration profiles in SFC requires using a model consisting of three sub-models: (i) a heat transfer model, (ii) a mass transfer model, and (iii) a model of the distribution of the mobile phase velocity.The model used here was successfully applied in modeling the temperature distribution inside the column and the elution of analytical concentration profiles during SFC [ 3 , 4 , 13 , 14 ].Here, the description is limited to the most important equations.A more detailed derivation of the isotherm model in the overloaded case is presented below.

The heat balance equation
The heat transport models for the column bed (1) and the column wall (2) are as follows: where C P,m , C P,s , and C P,w are the mobile phase, solid phase, and wall heat capacity, respectively (J m -3 K -1 ).The coefficient ε t is total porosity, u x and u r are the superficial mobile phase velocities in the x and r directions, respectively (m s -1 ), α is the mobile phase thermal expansion coefficient (K -1 ), λ eff is the effective bed conductivity (W m -1 K -1 ), λ w is the wall conductivity (W m -1 K -1 ), P is the local pressure (Pa), and T is the local temperature (K).
Using the above heat transport model, temperature profiles have been calculated for every applied flow rate and outlet pressure.The coefficient of thermal expansion of the mobile phase was calculated using the following equation: where ρ is the mobile phase density (kg m -3 ).Here, the effective thermal conductivity coefficient was calculated using the Zarichnyak correlations [15] : where the porosity, ε s , is the ratio of the volume of the solid phase in the bed to the geometric volume of the column, and λ s and λ elu are the solid phase and eluent conductivities, respectively (W m -1 K -1 ).

Mass balance equation
The equilibrium dispersive (ED) model is used to describe the mass balance over the column.The ED model in two dimensions (i.e., the axial and radial directions) can be written as:

∂ ∂r
r D ar ∂c ∂r (5) where c and q are the analyte concentrations in the mobile and stationary phases, respectively (g L -1 ), D ax and D ar are the apparent dispersive coefficients in the axial and radial directions, respectively (m 2 s -1 ), and u is a vector of the mobile phase velocity (see Eq. ( 13)).
The axial apparent dispersive coefficient, D ax , was calculated from the following equation [ 13 , 16 ]: where and where D L is the local dispersion coefficient (m 2 s -1 ), D m is the molecular diffusion coefficient (m 2 s -1 ), k ext is the external mass transfer coefficient (m s -1 ), d p is the particle diameter (m), and τ is the tortuosity coefficient.
The apparent radial dispersion coefficient, D ar , was calculated based on the plate height equation derived by Knox [17] : The external mass transfer coefficient, k ext , was calculated from the Wilson and Geankoplis correlation [18] : where and Sc = η ρD m (10) where Sh, Re, and Sc are the Sherwood, Reynolds, and Schmidt numbers, respectively.The dispersion coefficient, D L , was calculated from the following relationship [19] : (11) where γ 1 and γ 2 are the geometrical constants.It was assumed that γ 1 = 0.7 whereas γ 2 was estimated.The molecular diffusion coefficient, D m , was estimated from the Wilke-Chang equation [20] for the physicochemical conditions at the column inlet and outlet, and the average values of D m were taken for the calculations.The density and viscosity of carbon dioxide were taken from REFPROP 10 [21] .
The tortuosity parameter, τ , was calculated from the following correlation [22] :

Mobile phase velocity distribution and pressure distribution
The local value of the mobile phase velocity was calculated according to the following equation [23] : (13) where ( ρ/ η x ) denotes the average value of the ρ/ η ratio at a given axial position: and u o and ρ o are the mobile phase superficial velocity and the density, respectively, calculated at the column inlet.The pressure along the column was calculated using the correlation derived from the Blake, Kozeny, and Carman correlation: where ξ is an empirical parameter considered to equal 150, whereas the pressure drop was adjusted by estimating the external porosity.
The radial mobile phase velocity was calculated from the continuity equation:

Adsorption isotherm model
In this work, we applied the Langmuir model based on an equation proposed by Martire and Boehm [24] .According to these authors, the dependence of the retention factor on the reduced density and reduced temperature can be written as follows: where ρ R and T R are the reduced density and temperature, respectively.Parameters c 0 to c 4 are adjustable parameters to be estimated based on experimental data.Taking into account the definition of the retention factor, the linear isotherm model can be written as: The above equation was successfully used in predicting retention and concentration profiles at analytical scale [ 3 , 4 , 14 ].Including the first-order dependency of the reduced temperature and density on the retention factor, the following equation for nonlinear chromatography with Langmuir adsorption can be written as: where q s is the saturation capacity.

Instrumentation and procedure
The SFC system was an ACQUITY UPC2 system (Waters Corporation, Milford, MA, USA) equipped with a PDA detector, an autosampler with 100-μL loop was used for all injections.The Waters system utilizes an injection method called the "Partial Loop with Needle Overfill" which also inject a small amount of weak needle wash , used to transfer the sample liquid from a vial to the injector loop.In this case, MeOH was used as weak needle wash .The mass flow of the eluent was measured using a CORI-FLOW M12 Coriolis mass flow meter (Bronkhorst High-Tech B.V., Ruurlo, Netherlands) after the mixer.The pressures at the inlet and outlet of the column were measured using two EJX530A absolute pressure transmitters (Yokogawa Electric Corporation, Tokyo, Japan) connected to the column using a tee.The temperature of the column wall was measured at 20, 50, and 80% of the column length using three PT-100 four-wire resistance temperature detectors with an accuracy of ±0.2 °C (Pentronic AB, Gunnebo, Sweden).The temperature and pressure signals were logged using a PT-104 data logger from Pico Technology (St.Neots, UK).
All experiments were conducted at a set column temperature of 25 °C (298.15K).Three set instrumental outlet pressures were considered: 90, 110, and 150 bar.The flow rate of the mobile phase at all used outlet pressures ranged from 0.5 to 2.5 mL min -1 with 0.5 mL min -1 interval increments.The highest system pressure did not exceed 400 bar.In this condition, the chromatographic system was working in a subcritical state of carbon dioxide.The column with small-diameter particles and a variable mobile phase flow rate allowed us to obtain the experimental data with both small and large significant changes in the properties of the mobile phase inside the column.In this context it should be mentioned that the mobile phase used in this study as neat carbon dioxide with low temperature as well the use of the column with small-diameter particles are not typical for preparative SFC.However, it fully enabled the achievement of the main aim of this work which is the numerical modeling of overloaded elution concentration profiles by considering the changes in mobile phase properties in the axial and radial directions.
The solute was dissolved in toluene at concentrations of 0.75 g L -1 and 20 g L -1 for analytical and preparative injections, respectively.The injection volume was 2 μL in the analytical mode.Overloaded elution profiles were investigated by injecting 5, 10, and 20 μL of phenol.Competition for the active site between the diluent and the solute is limited and it is not included in the simula-tions.Two replicates of the injections were always made.The column was equilibrated for at least 40 min after changing the mobile phase flow rate to achieve steady-state conditions.The analytical and overloaded profiles were detected at 210 and 230 nm, respectively.The signal of the eluted overloaded profiles was converted to concentration using calibration curves determined separately for all considered mobile phase flow rates.The calibration curves were derived based on three peaks and for the conditions prevailing in the autosampler and at the outlet of the column.Additionally, the concentration profiles were adjusted to maintain the mass balance.

Method of estimating model parameters and simulations
To accurately calculate concentration profiles given the isenthalpic expansion of the mobile phase, first, the heat and pressure profiles inside the column must be calculated (see section 2.1 ).The unknown parameters of the model (i.e., external heat transfer coefficient, external porosity, and solid phase conductivity) were estimated using the column inlet pressure and the column wall temperature measured 8 cm from the inlet at a flow rate of 2.5 mL min -1 and set outlet pressure of 110 bar.Under these conditions, we expect the greatest changes in the mobile phase properties.
The adsorption isotherm model parameters were estimated in two steps.First, the analytical version of the isotherm model ( Eq. 19 ) was used.The model parameters were estimated from the phenol retention times under three different conditions, as follows: (i) a set outlet pressure of 90 bar and a flow rate of 0.5 mL min -1 , (ii) a set outlet pressure of 110 bar and a flow rate of 2.5 mL min -1 , and (iii) a set outlet pressure of 150 bar and a flow rate of 1.5 mL min -1 .Also, the true value of the geometric constant in the correlation of the dispersion (see Eq. 11 ) had to be determined.This was done based on the analytical peak recorded at a set outlet pressure of 150 bar and a mobile phase flow rate of 1 mL min -1 .Second, the saturation capacity was estimated from one overloaded elution concentration profile obtained at a set outlet pressure of 110 bar and a mobile phase flow rate of 2.5 mL min -1 .When estimating the saturation capacity, the parameters determined in the first step were kept constant.
When calculating the concentration profiles, first the axial and radial temperature gradients in steady state were calculated as a function of the physicochemical properties of the mobile phase and its flow rate.Then the mass transfer model was solved using the temperature and pressure obtained from solving the heat transfer model.The physicochemical properties of the mobile phase, at each stage of integrating the model equations, were calculated using REFPROP v. 10 [21] .To discretize the spatial derivatives of the model, the orthogonal collocation on finite elements method was used [26] .To solve the system of ordinary differential equations, the Adams-Moulton method implemented in the CVODE procedure was used [27] .In all estimations of the unknown parameters, the inverse method with the Levenberg-Marquardt algorithm was applied [28] .

Results and discussion
The aim of this study is to numerically model overloaded elution concentration profiles in SFC by considering the changes in mobile phase properties in the axial and radial directions.This was done in two steps, which are presented below.First, in section 4.1 , the heat transfer model is calibrated and validated.Moreover, the influence of the temperature and pressure distribution on the properties of the mobile phase is discussed.In section 4.2 , the results of calculating the overloaded concentration profiles using the temperature and pressure distributions are presented and discussed.To accurately model elution profiles in SFC under non-isopycnic conditions, the expansion of the mobile phase also needs to be considered.The viscous friction will heat the column, whereas the expansion of the mobile phase will cool the column.Therefore, the heat and pressure profiles inside the column have to be calculated.For this purpose, the two-dimensional heat transfer model was used (see section 2.1 ).Several unknown parameters must be determined: the external porosity, the effective thermal conductivity of the solid phase, and the external heat transfer coefficient (see section 3.3 for more details).In Table 1 , these unknown parameters are presented using typical values found in chromatographic systems.The external heat transfer coefficient is used in boundary conditions to solve Eq. ( 1) and capture the heat transfer between the column oven and the column wall.Here, the external heat transfer coefficient is 22.85 W m -2 K -1 for the air-controlled temperature mode.The external porosity is 0.3462 based on the Blake, Kozeny, and Carman correlation and the permeability coefficient is 150.The effective conductivity of the solid phase was estimated to 0.7998 W m -1 K -1 , see section 3.3 for more details.The matrix of the particle is inorganic silica with a conductivity of 1.4 W m -1 K -1 .However, as observed here the ligands of the stationary phase are 2-(3-propoxy) ethane-1,2-diol silanes, have significantly reducing the heat conductivity of the stationary phase.
The heat balance and estimated velocity distribution were validated by comparing the calculated pressure at the column inlet and the calculated wall temperature in three positions (i.e., 2, 5, and 8 cm from the column inlet) with the corresponding experimental values.The calculated and measured temperatures were in good agreement in all runs.The maximum relative error was less than 0.13%, but the error did not exceed 0.05% in most cases.This indicates that the maximum difference was less than 0.41 K, but the difference did not exceed 0.2 K in many cases.The maximum relative error of the calculated inlet pressure was 1.84%.The maximum difference in the calculated pressure was 2.76 bar, whereas the difference in over 50% of runs was about 1 bar or even significantly less (for more details, see Tables S.1 and S.2 in Supplementary material).This good agreement confirms the validity of using the heat transfer model and the correctness of the pressure and temperature distribution calculations in the studied case.
As expected, the greatest changes in the properties of the mobile phase occur at the lowest considered backpressure and the highest mobile phase flow rate.Fig. 1 a presents the calculated temperature of the column bed and wall at a mobile phase flow rate of 2.5 mL min -1 and a set outlet pressure of 90 bar.The temperature drop along the column exceeds 5 K while the temperature difference between the column center and the inner column wall is higher than 1 K close to the column outlet and reaches the highest value equals 1.70 K at the end of the column.At a set outlet pressure of 150 bar, the temperature drop is smaller: we observe 2.54 K and 0.91 K for the temperature drop along the column and the temperature difference in the radial direction, respectively (see Fig. 1 b).One may be surprised about the small observed radial temperature gradient in these cases.However, it can be explained by comparing the column's geometry and heat transport.Here, the ratio of column length ( L ) to its radius ( R ) is equal to about 60.Let us assume that the heat Q transported in the axial ( Q x ) and radial direction ( Q r ) are equal.We know that R is equal to about 120, so in other words, the changes in temperature in the radial direction as compared to axial direction would be 120 times smaller.
Fig. 2 a shows the density of the mobile phase along the column at a mobile phase flow rate of 2.5 mL min -1 and a set outlet pressure of 90 bar (note that the column wall is not plotted as in Fig. 1 ).Here, we observe a density drop of 116 kg m -3 along the column and of 10 kg m -3 in the radial direction at the column outlet.Like as the temperature changes, at a higher set outlet pressure of 150 bar, the density changes inside the column are smaller, being 98 kg m -3 and 4 kg m -3 in the axial and radial directions, respectively (see Fig. 2 b).Note that the scale differs between Fig. 2 a and 2 b due to the large difference in the density at the 90 and 150 bar set outlet pressures.It should also be noted that the calculated pressure drop (not presented in the figure) at the 2.5 mL min -1 volumetric flow rate exceeds 200 bar and increases with increasing outlet pressure, i.e., 203 and 221 bar at the 90 and 150 bar set outlet pressures, respectively.In this region, CO 2 is a compressible liquid (i.e., subcritical condition).Here, the viscosity increases as the pressure increases and decreases as the temperature increases.Compare Fig. 2 in which the viscosity is plotted at set outlet pressures of (c) 90 and (d) 150 bar, respectively.The mobile phase properties in this region are thus similar to those of typical liquids.
From Eq. ( 13) follows that the properties of the mobile phase influence the mobile phase velocity, which is no longer constant across and along the column and can have a significant influence on the eluted profiles.Fig. 3 shows the calculated superficial velocity inside the column at a 2.5 mL min -1 volumetric flow rate and set outlet pressures of (a) 90 and (b) 150 bar, respectively.The thermal expansion of the mobile phase results in the increase of the mobile phase velocity and thus the increase of the concentration profile movement.The increase in velocity is closely related to the density drop, see Fig. 2 a and 2 b.In the cases studied, the superficial velocity increases from the inlet to the outlet of the column by 12 and 10% at set outlet pressures of 90 and 150 bar, respectively.The radial gradient of the velocity due to the temperature gradient is also observed, see Fig. 3 a and 3 b.The velocity increases from the colder center of the column to the warmer wall region, except in the region near the column inlet where the opposite tendency is observed.

Overloaded concentration profiles
Above, the temperature and pressure distributions, as well as the velocity, viscosity, and density of the mobile phase, were calculated; here they are used to calculate overloaded concentration profiles for solving the two-dimensional mass transport model.One of the most important things in simulating concentration profiles is to use an appropriate adsorption isotherm model -in this case, the modified Langmuir isotherm model, which also considers the reduced density and reduced temperature (see Eq. [19] ).The model parameters were estimated in two steps (see section 3.3 ) and are presented in Table 2 .It is worth noting that the first three parameters, i.e., c 1 to c 3 in Eq. ( 17) , are enough to forecast the retention times in the considered experimental condition.It follows from the calculated confidence level that the most important term in the isotherm model is parameter c 2 , which is connected with the reduced density (see Table 2 ).In other words, the density of the mobile phase is the major factor that controls solute retention.Moreover, the values of parameters c 1 and c 2 are positive and negative, respectively.This indicates that the solute retention decreases at a constant temperature with increasing density and that the retention increases with decreasing temperature at constant density; both these observations are typical of SFC.In Fig. 4 , the retention factor for phenol is plotted as a function of the density and temperature of the mobile phase.
The γ 2 parameter in Eq. ( 11) had to be determined.This parameter was estimated based on the analytical peak obtained at a mobile phase flow rate of 1 mL min -1 and a set outlet pressure  of 150 bar.This peak was chosen because the smallest temperature and density gradients are observed in this condition.The estimated value of the parameter was validated based on other analytical peaks and good agreement was indicated.The final value of γ 2 was 17.8, which would seem to be very high.However, the analytical peaks in all considered runs are not sharp.It should also be emphasized that this column is not a common column available for sale.It was specially ordered and the particles are smaller in diameter than those used in common SFC columns.The increased dispersion could be because of imperfectly spherical particles and/or heterogeneous packing of the column.Now only the saturation capacity in the adsorption model needs to be determined (see section 3.3 for details).The saturation capacity was determined to be 252.5 g L -1 (see Table 2 ).
Finally, the estimated adsorption isotherm model was used to predict band profiles for all other experimental conditions than those used during model calibration.Fig. 5 shows the experimental and simulated elution profiles for 5-, 10-, and 20-μL injections of the sample with a concentration of 20 g L -1 and at mobile phase flow rates of (a) 0.5 mL min -1 , (b) 1 mL min -1 , (c) 2 mL min -1 , and (d) 2.5 mL min -1 .The set outlet pressure is 110 bar (for comparison with 90 and 150 bar see Supplementary Material Figs.S1 and S2, respectively).The agreement between the calculated and experimental concentration profiles is good.In all considered experiments, the proposed model can forecast the overloaded concentration profiles well.It can be observed that the shock layers of the peaks are not as sharp as predicted.The front dispersion increases with decreasing injection volume and increases slightly with increasing mobile phase flow rate.To investigate whether this effect is due to dispersion, elution concentration profiles were calculated for different dispersion coefficients, see Fig. 6 .The calculations clearly show that the distortion was caused by the dispersion.
As stated above, good agreement between the calculated and experimental concentration profiles was obtained.Only a small shift in retention times was observed, caused by imperfect prediction of the retention times using a linear isotherm model.However, including the two rejected parts of linear isotherm model parameters did not improve the model's ability to predict the solute retention.However, the predicted heights and widths of the calculated elution profiles are in very good agreement with corresponding experimental profiles, even though the saturation capacity was estimated based on one peak.This indicates the correct calculation of the saturation capacity and, especially, of the association equilibrium constant, which changes along the column.As can be seen, the density drops along the column exceed 100 kg m -3 at Fig. 5. Comparison of experimental (dotted line) and simulated (solid line) concentration profiles at mobile phase flow rates of (a) 0.5 mL min -1 , (b) 1 mL min -1 , (c) 2 mL min -1 , and (d) 2.5 mL min -1 .The set pressure at the outlet of the column is 110 bar.Phenol was used as a solute with a sample concentration of 20 g L -1 .The injection volumes were 20, 10, and 5 μL. .The experimental profile is plotted with a blue dotted line while the calculated profiles are plotted with a green dash dot line (for γ 2 = 1), with a red dash line ( γ 2 = 10) and a black solid line ( γ 2 = 17.8), respec- tively.The set outlet pressure is 90 bar while the mobile flow rate is 2.5 mL min −1 .The injection volume is 5 μL.(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)Fig. 7. Calculated solute concentration profile inside the column just before the zone elutes at a mobile phase flow rate of 2.5 mL min -1 and a set outlet pressure of 90 bar.Otherwise, the conditions are as in Fig. 6 d.
the highest mobile phase flow rate, causing a drastic change in the retention factor along the column.With such a density drop, the retention factor can increase by approximately 80%.Without taking into account the adsorption isotherm calculations at each point along the column, the peak retention and its shape could not be correctly simulated.
As shown in this study, two-dimensional column models allow a more detailed investigation of the movement of elution profiles, since they account for both the axial and radial directions.This is essential because the temperature, viscosity, and mobile phase velocity vary in both the axial and radial directions.Fig. 7 shows the calculated elution concentration profile near the column outlet at the highest mobile phase flow rate and a set outlet pressure of 90 bar.In Fig. 7 one does not observe any difference in solute migration velocity in the radial direction, even at the lowest outlet pressure, at which the largest density drops occur.The underlying reason for this is that the effects of radial temperature and density gradients and mobile phase velocity gradients cancel each other out.However, it is worth emphasizing the phenomena occurring in SFC even though they do not have significant effects in this particular case.Let us discuss this in more detail for the second half of the column -that is, the part from the middle of the column to the outlet in the axial direction.The temperature increases from the center of the column to its wall (see Fig. 1 ), while the density of the mobile phase is more or less inversely proportional to the temperature (see Fig. 2 ).Thus, the density decreases from the column center to the column wall.The velocity of the mobile phase is higher close to the wall and not in the center of the column, as is generally observed in LC systems (see Fig. 3 ), so the velocity gradient will make the solute move faster in the wall region.However, the speed of the solute migration also depends on the retention factor, which depends on both the temperature and density of the mobile phase (see Fig. 4 ).From Fig. 4 we can observe that, in this particular case, the density was the most important factor, and that temperature generally only slightly affected the retention.As stated above, the density of the mobile phase is highest in the center of the column, so the retention factor is the lowest in the center of the column.In the first half of the column, from the inlet to the middle, the radial gradients discussed above are smaller.For the temperature, the gradient is opposite to that at the end of the column, because the temperature of the column wall is lower than the temperature of the mobile phase entering the column (see Figs.

Comparison of one and two-dimensional models
In this section, we compare the solution of the proposed twodimensional column model (2D) with its one-dimensional version (1D).The advantage of a 1D model is that the calculations could be performed several dozen times faster as compared with the 2Dmodel, which is of great importance e.g. if model parameters are estimated with the inverse method or process is optimized.
In the 1D model we make the following assumptions: 1) The column temperature is equal to oven temperature, -we neglect the drop of temperature in the axial direction, which seems to be justified in this case because temperature drop was generally less than 5 °C 2) Gradients of physicochemical parameters do not exist in the radial direction 3) The pressure inside the column is decreasing linearly and only in the axial direction, from column inlet to the column outlet -this assumption is justified by the observation made among others in [29] 4) Mass flow in any cross-section of the column is constant and equal mass flow at the inlet to the column.
We can now write the mass balance equation in 1D-model as: Having pressure distribution we calculate density as a function of position along a column on the base of REFPROP 10.
According to assumption, 4 flow rate is calculated from the equation: The axial dispersion coefficient, D ax , is calculated as before.Fig. 8 presented the comparison of peaks calculated at different flow rates (0.5 mL min −1 and 2.5 mL min −1 ) and back pressures (90 bar and 110 bar) with both models.The red dash line represents the solution with the 1D model and the black line with the 2D model.As can be seen, the peaks are shifted slightly -the difference does not exceed 2%.The difference in the position of calculated concentration distribution depends on changes in flow rate and temperature along the column.Larges deviations were observed at the highest flow rate, see Fig. 8 b and 8 d as compared to the lower flow rate, see Fig. 8 a and 8 c.Because the flow rate is increasing and the temperature is decreasing their influence on retention is reciprocal.However, in the 1D model, the temperature changes are ignored, which can result in observed changes in peaks positions.It should be also noticed that axial velocity, calculated with Eq. ( 21) is not the same as calculated in the 2D-model, see Eq. (13) .
The differences between elution zones calculated with 1D and 2D models will depend on the temperature and pressure dependency of adsorption.The aim of this study is not to formulate a rule of thumb under which the 2D model could be replaced with a simple 1D model, without serious loss of accuracy.Therefore, the simpler 1D model should be used with the great caution.

Conclusion
This study shows that it is difficult, even impossible, to obtain a complete census of all interactions and thus to fully understand how temperature, density, and velocity gradients affect separation in SFC, without the aid of numerical calculations.
The overloaded concentration profiles in SFC were predicted using combined two-dimensional heat and mass transfer models connected with the model that allowed the calculation of the superficial velocity of the mobile phase.The isotherm model was derived based on a relation used for predicting the retention factor in SFC.The equilibrium constant in the isotherm model was a function of the temperature and density of the mobile phase.Thus, the influences of mobile phase properties on retention and of the movement of the elution profiles due to velocity distribution on the shape of the concentration profiles and elution times were taken precisely into account.
The prediction was made using the experimental overloaded band profiles of phenol eluted with liquid carbon dioxide on a col-umn packed with 2.2-μm-diameter silica gel particles with diol ligands.The experimental data consisted of elution profiles recorded at three injection volumes, five mobile phase flow rates, and three set outlet pressures.Moreover, the mass flow, the pressure at the column inlet and outlet, and column wall temperature were controlled.The parameters of the heat transfer model, as well as the isotherm model parameters were estimated based on experimental data using the inverse method.The calculated temperature distribution and pressure along the column agreed very well with the corresponding experimental values.
Depending on the conditions of the chromatography process, insignificant to large changes in the mobile phase temperature and density were observed.At the highest mobile phase flow rate (2.5 mL min -1 ) and lowest set outlet pressure (90 bar), the density drop along the column was 116 kg m -3 while the density changes in the radial direction at a distance of 1.5 mm reached 11 kg m -3 ; these correspond to the temperature drops, which were 5 and 1.7 °C in the axial and radial directions, respectively.The mobile phase temperature and density drop resulted in other changes in mobile phase properties.As a result of these changes, the superficial velocity increased from the column inlet to outlet by 12% at the highest considered volumetric flow rate.The retention factor in such conditions increased by up to 80%.Increasing the set outlet pressure slightly decreased the radial and axial gradients, while reducing the mobile phase flow rate to 0.5 mL min -1 caused the column to operate under more or less flat profiles of temperature, density, and mobile phase velocity distribution.
Concentration profiles calculated using the two-dimensional ED model taking into account the changes in mobile phase properties and velocity distribution agreed very well with the corresponding experimental concentration profiles.A small shift in the elution times was observed.However, the shape of the concentration profiles was perfectly maintained in all considered runs.Moreover, distortion of the shock concentration was observed as a result of the increased dispersion.The radial gradient was not observed to influence the concentration profiles.
We showed that the applied methodology and mathematical model, including the isotherm model derived for overloaded SFC, lead to good agreement between theoretical and experimental results.The model allowed us to predict overloaded concentration profiles with good accuracy over a wide range of mobile phase changes.In the end, we also formulated a simpler onedimensional column model which with caution in some cases can be used for modeling overloaded SFC.We must stress that the twodimensional column model is general and will work in all cases and is therefore used in this study.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Temperature distribution (in Kelvin) calculated at set pressures at the column outlet: (a) 90 bar, (b) 150 bar.The mobile phase flow rate is 2.5 mL min -1 .Radius 0 is the center of the column.The dotted line represents the inner column wall.

Fig. 4 .
Fig. 4. Calculated retention factor of the solute as a function of density at four different mobile phase temperatures (19, 21, 23, and 25 °C).The retention factor was calculated based on the estimated isotherm model parameters.

Fig. 6 .
Fig. 6.The simulated concentration profile calculated for different values of the parameter gamma 2 ( γ 2 , in Eq. 11 ) .The experimental profile is plotted with a blue 1 a and 1 b).

Table 1
Physicochemical properties of the column included in the study.

Table 2
Estimated parameters of the isotherm model.