Adsorption Isotherm Modelling of Water on Nano-Tailored Mesoporous Silica Based on Distribution Function

Abstract

A new model of adsorption isotherms Type IV and V is proposed as a basis for theoretical calculations and modelling of adsorption systems such as adsorption heat storage and heat pumps. As the current models have decent yet limited applicability, in this work, we present a new combined model with universal use for micro-mesoporous silica/water adsorption systems. Experimental measurement of adsorption isotherm of water onto seven different samples of micro and mesoporous silica and aluminium-silica were used to fit new adsorption models based on a combination of classical theories and a distribution function related to the pore-size distribution of the selected materials. The fitting was conducted through a repeated non-linear regression using Trust Region Reflective algorithm with weighting factors to compensate for the scalability of the adsorption amount at low relative pressure with optimization of the absolute average deviation fitting parameter. The results display a significant improvement for most of the samples and fitting indicators compared to more common models from the literature with average absolute deviation as low as AAD = $0.0025$ $\mathrm{g}$ $\mathrm{g}$⁻¹ for material with maximum uptake of q = $0.38$ $\mathrm{g}$ $\mathrm{g}$⁻¹. The newly suggested model, which is based on a combination of BET theory and adjusted normal distribution function, proved to bring a higher degree of precision and universality for mesoporous silica materials with different levels of hydrophilicity.

1. Introduction

Heat driven thermodynamic cycle is as an alternative to the common mechanical vapour compression cycle, where an abundance of waste heat or low potential heat (low-temperature difference) is available. Adsorption cycles can work with minimum consumption of electricity provided sufficient heat source such as solar thermal [1] or co-generation [2,3,4] providing heat or cold from renewable sources or waste, based on the application. The adsorption phenomena works as a driving force for these systems, and adsorption isotherm models are an essential tool to describe and predict their behaviour. The adsorption models are readily used in simulations of adsorption heat storage (AHS) [5,6] and adsorption heat pump (AHP) [7,8,9,10,11] systems. The advantage of such an approach is a low necessary amount of the actual sample and faster evaluation of the key properties of the adsorption system.

Majority of the research of the adsorption heat-based technologies involves AHP thanks to their more efficient performance and significantly less mass of the deployed system as well as the lower amount for the expensive adsorption materials. However, both the AHP and the AHS systems can be modelled based on the isothermal data with good accuracy and analysis results of the key performance points. The most common models are then deployed based on the adsorption material characteristics and their practicality. Dubinin-Astakhov [8,9] and Tóth [5,11] adsorption models are usually used to describe materials with high adsorption potential such as silica gel, which generally shows high water affinity as Type I isotherm. Helaly et al. [5] used Tóth adsorption model for performance analysis of water adsorption on silica gel and zeolite 13X-based adsorption storage systems and achieved a good agreement with the actual experimental data. Likewise, Schmidt et al. [6] showed on the example of zeolite 13x/water adsorption pair that a proper system analysis of its performance can be obtained just from the adsorption model considering its physical validity. [11] used in their work uncommon Mahle’s isotherm model [12] to fit the data of a Type V isotherm. However, because the Mahle’s model does not properly consider the surface adsorption and focuses mainly on the capillary condensation part, it can be technically successfully deployable only on Type V isotherms. Hence, to jointly describe also the adsorption characteristics of silica gel, ferroaluminophosphates FAM Z01 and FAM Z02, and metal-organic frameworks MIL 101(Cr) and NH${}_2$-MIL 125(Ti) they used Tóth isotherm [13] next to the Mahle’s model. For more complex isotherm showing Type IV isotherm some unorthodox high polynomial functions have been documented as well [10]. The problematic of current isotherm models used to describe the adsorption uptake in the adsorption cooling applications were also identified by Mohammed et al. [14] with the introduction of some corrections to the most common adsorption material, silica gel. That is, however, mostly for microporous materials or materials with high water affinity showing uncomplicated isotherms curves of Type I.

Accordingly, the mathematical models describing the water adsorption isotherms on mesoporous materials are a subject of many research publications [15,16,17,18,19]. However, on the contrary to the classical models which exhibit simple parabolic or hyperbolic curvature of the adsorption isotherms based on the pure surface adsorption, porous materials display frequently complex adsorption behaviour which can be described by those models only partially. That means that these can provide us only with a limited precision or can be applied only to a specific region of the isotherm [20]. The sigmoidal (or S-shaped) adsorption isotherm, typical for the mesoporous materials, is usually identified by the International Union of Pure and Applied Chemistry (UIPAC) as Type IV or Type V, plus subtypes (a) and (b), depending on the immediate surface characteristics [21]. The sigmoidal shape shown in Figure 1 is typically caused by secondary adsorbate–adsorbate interactions such as capillary condensation in mesopores. Different shapes of the individual types are then given by their specific adsorbent–adsorbate relation characteristics. Type IV corresponds to materials with a surface more susceptible or accessible to adsorption of a given adsorbate, while Type V describes adsorption processes burdened by certain constraints. Those constraints can be of chemical characteristics, such as lower affinity towards the selected adsorbate (e.g., silica without hydroxyl groups with water), or of physical characteristics in case of a small opening window to the pores of some zeolite materials [22,23,24,25] resulting in higher activation energy. The Type IV isotherm is further divided to Type IV(a) with an adsorption/desorption hysteresis, and Type IV(b) without hysteresis [21]. A particular case represents the stepwise Type VI isotherm typical for some layer-by-layer adsorption on highly uniform nonporous surfaces [16]. However, we can also observe similar behaviour at significantly higher relative pressures for water adsorption on some mesoporous materials with a particular pore size distribution like the TMPS-2.7A reported in our previous findings [26].

While the simplest of this selection is possibly the Type V isotherm with only one rising edge in the small range of the relative pressure, the most complicated adsorption systems are described by Type VI. Although the adsorption systems of Type IV and Type V isotherms are well known and abundant, thanks to the nature of the adsorption process on mesoporous material, the theoretical models are still difficult to universally apply and new models are continuously developed. The modelling attempts of the isotherm adsorption curves on mesoporous materials based only on the simple empirical expressions and fundamental adsorption models such as Henry’s [27], Langmuir’s [28], Freundlich’s [29], Temkin’s [30], BET [31] and their algebraic combinations [15] usually fail due to the complicated shape of the actual isotherm or insufficient accuracy; or can be used only on selected examples of adsorbent–adsorbate pairs.

2. Sigmoidal Adsorption Isotherm Models

2.1. Classical Models: Development of Multilayer Adsorption Models

The BET isotherm (Brunauer, Emmett, Teller) was one of the first theories widely recognized for dealing with considerations of adsorbent–adsorbate and adsorbate-adsorbate interactions. Built on Langmuir’s model, the BET model adds a theory of the layered structure and interaction between the layers [31]. As a result that it follows in the low relative pressure of the Langmuir model, it also sets the maximum amount $q_m$ of the sorbate to cover the bare surface of the adsorbent completely. The BET isotherm is written in its simple form as:

$$\theta =\frac{q_e}{q_m}=\frac{C\frac{P}{P_0}}{\left(1-\frac{P}{P_0}\right)\left(1+\left(C-1\right)\frac{P}{P_0}\right)}$$

(1)

where $\theta$ is the relative adsorption amount, $q_e$ is the adsorption amount at the equilibrium, $q_m$ is the adsorption capacity of the surface layer (maximum adsorption in case of Langmuir’s model) and C is the BET constant describing the adsorbent–adsorbate interaction intensity. Thanks to a good fit with the experimental data of nitrogen adsorption at cryogenic temperatures for the majority of the materials, it is often used as a primary method to assess the surface area in its linear region [21,32], providing the validity $\theta =q_e/q_m$, $q_e<q_m$. If the adsorption follows the BET model, the resulting isotherm takes a linear form of Equation (2), and the BET constant C and the volume of adsorbate necessary for covering the surface $v_m$ can be obtained. Consequently, knowing the coverage factor (volume per area) of the nitrogen for the given material, the surface area can be calculated.

$$\frac{P}{q_e\left(P_0-P\right)}=\frac{1}{q_{m}C}+\frac{C-1}{q_{m}C}\frac{P}{P_0}$$

(2)

Although the BET provides well-fit parameters in the lower relative pressure, the higher relative pressures $P/P_0>0.35$ are burdened by a significant error [33]. The cause for this may vary; however, one of the explanation can be the underestimating of the adsorption enthalpy for second and higher layers in correlation to the change of the relative pressure in the model. For this reasons the original BET model was modified originally by Anderson [34], and later derived again through different approach by Brunauer et al. [33], by adding another constant to the relative pressure; $P/P_0\to K_a\left(P/P_0\right)$. This modification was also introduced and derived independently by de Boer [35] and Guggenheim [36], and therefore it became broadly known as the Gugenheim, Anderson, de Boer (GAB) isotherm [37,38]. By practically limiting the relative pressure factor reaching unity, $K_a<1$, the model is prevented from reaching infinite adsorption at full saturation $P=P_0$ and an apparent equilibrium can be obtained, nevertheless, leaving the final adsorption infinite. However, even with this extension, the BET model can be used only to define Type II and Type III isotherms when the system is reaching full saturation. For more complex systems the modified BET adsorption model fitting follows only the first phases of the adsorption process.

The adsorption of multiple adsorbate molecules or atoms can also be simply described extending the previous formula for monolayer adsorption as introduced by Klotz et al. [39]:

$$\begin{array}{cc}\hfill S+A& \leftrightarrow SA\hfill \\ \hfill SA+A& \leftrightarrow S{A}_2\hfill \\ \hfill S{A}_2+A& \leftrightarrow S{A}_3\hfill \\ \hfill \cdots \end{array}$$

(3)

Equation (3) implies that one adsorption site S can be a host for several adsorbate particles A until a saturation point or a maximum number of bonds is reached. The following equilibrium constants k are then naturally given as the ratio between the occupied and unoccupied sites together with free adsorbate particles:

$$k_n=\frac{\left[S{A}_n\right]}{\left[S{A}_{n-1}\right]\left[A\right]}$$

(4)

where n is the number of bounds (layers) that can be associated with one specific adsorption site as $n\le 1$. The adsorption fraction $\theta$ at any point is then given as:

$$\theta =\frac{k_1\left[A\right]+2k_1{k}_2{\left[A\right]}^2+\cdots +n\left(k_1{k}_2\cdots k_n\right){\left[A\right]}^n}{1+k_1\left[A\right]+k_1{k}_2{\left[A\right]}^2+\cdots +\left(k_1{k}_2\cdots k_n\right){\left[A\right]}^n}$$

(5)

Even though that Klotz et al. [39] used this model for ligands adsorption on proteins, this common approach was and is readily used by others for vapour adsorption as well. Hence, the adsorbate concentration $\left[A\right]$ in Equation (5) can be expressed as the relative pressure $P/P_0$. In a system with multiple layers, Equation (5) can become unnecessarily complicated and an approximation is usually considered. In the case of multilayer/multi-bond adsorption, we can simplify Equation (5) by assuming that the higher layers’ sorbate–sorbate interaction is constant; e.g., water adsorption on non-polar materials such as microporous activated carbon; and reduce the model to:

$$\theta =\frac{{\sum }_{n=1}^{i}n{k}_1{k}^{n-1}{\left(\frac{P}{P_0}\right)}^n}{1+{\sum }_{n=1}^{i}n{k}_1{k}^{n-1}{\left(\frac{P}{P_0}\right)}^n}$$

(6)

where, following Equation (5), $k_1$ is the constant describing the adsorption on the surface and k is the contribution from subsequent adsorptions connected to the original primary adsorption site S. Finally, the n factor is the association number corresponding again to the multilayer/multi-bond structure. The theory behind the original BET model [31] is very similar; however, it considers the properties of the second and above layers the same as the liquid state of the adsorbate, which practically leaves $k=1$. The derivation then naturally takes the path for the classical BET model which gives two results depending on considering the infinite number of layers (Equation (2)) or some limitations of the adsorption system; e.g., porous structure. In the latter case, the association number n is finite, or relatively small, and the BET equation will take a different form from the classical interpretation to accommodate for the space limitations:

$$\theta =\frac{Cx}{1-x}\frac{1-\left(n+1\right)x^n+n{x}^{n+1}}{1+\left(C-1\right)x-C{x}^{n+1}}$$

(7)

where C is again the standard BET constant. The factor x in Equation (7) in the original derivation is the typical relative pressure $P/P_0$. Many other isotherm models used a similar approach as well to the problematic of adsorption according to Equation (5) and reach the same or very similar conclusion after deriving the final formula. Similarly to GAB isotherm, an introduction of an additional constant was performed by several investigators to compensate for the surface-independent interactions giving:

$$x=K_a\left(P/P_0\right)$$

(8)

This interpretation was also achieved by Buttersack [19] by deriving Klotz’s [39] (Equation (6)), assuming a reduced influence of the nucleation points of the adsorption allowing:

$$\theta =\frac{{\sum }_{n=1}^{i}n{x}^n}{1+{\sum }_{n=1}^i{x}^n}$$

(9)

The $K_a$ from substituting Equation (8) into Equation (9) also has similar meaning with the previous expressions and in the case of vapour adsorption in a porous matrix corresponds to capillary condensation. That is, incidentally, also the case of identical expression of the $\zeta$-isotherm derived by Ward and Wu [40] in their attempt to develop a method for determining surface tension based on clustered adsorption. Buttersack [19] follows this approach as well; however, he suggests a semi-empirical linear extension to the Klotz model to better accommodate the condensation/sorption occurring at higher relative pressures.

$$\theta =\frac{Cx}{1-x}\frac{1-\left(n+1\right)x^n+n{x}^{n+1}}{1+\left(C-1\right)x-C{x}^{n+1}}+K_B\frac{P}{P_0}$$

(10)

where the C, n has the same meaning as in Equation (7) and x is substituted according to Equation (8), hence employing the apparent relative pressure, rather than the relative pressure itself. The newly introduced constant $K_B$ is obtained by non-linear regression of the experimental data. Naturally, this expression inevitably changes the original meaning of the constants from Equation (7), with the Equation (8) substitution, to a certain degree. Hence, even though, the results of the fitting show some enhancement of the statical criteria compared to the “plain” Klotz’s equation [19], the author himself encourages careful interpretation of the data together with high accuracy of the underlying data.

2.2. Cooperative Multi-Molecular Sorption and Rutherford’s Model

As a result that the standard models are somehow insufficient in the description of the S-shaped isotherms, attempts with the compliance of the general models are typically pursued on a case basis with appropriate adjustments to the adsorption theories. Based on the classical Langmuir concept, Malakhov and Volkov developed the Cooperative Multi-molecular Sorption (CMMS) model [41] mostly to describe the progress of the adsorption close to the Type V isotherm. In their work, they focus on adsorption of simple alcohols (MeOH, EtOH, iso-PrOH) on a microporous polymeric structure, which can be described as a network of diffusion channels or nano-voids with a diameter of roughly 1.4 nm [41]. According to the CMMS theory, the adsorption progress alongside adsorption channels, the polymer chains and the adsorbed molecules promote additional adsorption mainly in the adjacent empty and open sites, thus, cooperation. The following equation of the adsorption fraction $\theta$ then describes the process based on rate (filling) constants ($K_{as},K_0,K_1$), each corresponding to the effective adsorption on the unoccupied sites in the proximity of the primary sorption sites:

$$\theta =\frac{q}{q_0}=\frac{K_0\frac{P}{P_0}}{\left(1-K_{as}\frac{P}{P_0}\right)\left(K_0\frac{P}{P_0}+W^2\left(1-K_{as}\frac{P}{P_0}\right)\right)}$$

(11)

The Equation (11) is given for vapour adsorption and the adsorption is expressed as a fraction of the maximum adsorption $q_0$ and adsorption q at the given pressure P through the relative pressure $\frac{P}{P_0}$ at the equilibrium point. An important factor is the interaction parameter W employing the Ising model, which explains the cooperative effect of the occupied sites. The parameter W is given as:

$$W=\frac{1}{2}\left[1-K_1\frac{P}{P_0}+\sqrt{{\left(1-K_1\frac{P}{P_0}\right)}^2+4K_0\frac{P}{P_0}}\right]$$

(12)

where $K_0$ and $K_1$ are the filling constants for adjacent and adjacent-adjacent primary sites, respectively, in the direction of the adsorption chain. The adsorption rate $K_{as}$ from Equation (11) is then intended for the adsorption rate in the lateral direction (increasing number of layers). Furthermore, if we consider that the cooperative effect is zero or negligible ($K_0/K_1=1$), the value of the interaction parameter becomes $W=1$ and Equation (11) is transformed into the extended BET isotherm as described by Anderson [34] and Brunauer et al. [33]. Then withal, if the lateral adsorption rate is $K_{as}=1$, Equation (11) becomes the standard BET [31].

Rutherford [42] realized that the CMMS could be used on some other porous materials as well provided the similarity of the systems based on the strong sorbate–sorbate interactions and weak sorbent–sorbate affinity, also observed for Type V isotherms. Rutherford [42] used the CMMS on an example of adsorption pair of microporous carbon/water. In the following work, Rutherford [43] reasons that the adsorption of water onto microporous carbonous surfaces consists of two modes; one which considers water bonding to the functional groups, and the second when water adsorbs within the micropores. As a result that the original CMMS cannot comprehend such duality, Rutherford [43] proposed an extended form using two CMMS equations besides. The extended Rutherford’s adsorption model was then also successfully applied in other works using microporous carbon [44]. The combined form can be formally written as:

$$q=q_{satL}CMM{S}_1+q_{sat}CMM{S}_2$$

(13)

where q is the total adsorbed amount from both modes, $CMM{S}_1$ and $CMM{S}_2$ are the respective forms/contributions of the adsorption equations in the individual modes, and $q_{satL}$ and $q_{sat}$ are the adsorption capacities for respective modes. Upon considering the adequate values of the rate factors $K_{as},K_0,K_1$, the final form can be written:

$$q=\frac{q_{satL}{b}_L\frac{P}{P_0}}{1+b_L\frac{P}{P_0}}+\frac{q_{sat}{K}_0\frac{P}{P_0}}{K_0\frac{P}{P_0}+W^2}$$

(14)

The first part of Equation (14) represents the first mode; the interaction of the adsorbate with the functional groups on the surface, and the second part describes the adsorbate molecules adsorption in micropores; i.e., the second mode. The considerable simplification of the first mode is allowed by the fact that the adsorption to the functional groups is direct and occurs only in one layer. Therefore, the lateral adsorption is not considered and the adsorption rate $K_{as}=0$. Furthermore, considering that the adsorption is dependent on the positions of the functional group and adsorption on the adjacent positions has the same probability, the $K_0=K_1$ and, thus, the $W=1$ and the CMMS takes formally the shape of the Langmuir isotherm (L index). The constant $b_L$ from Equation (14) is, therefore, the standard Langmuir affinity constant. For the second mode, the adsorbate–adsorbate interactions have to be considered again as well, giving $K_0\ne K_1$, while the lateral adsorption is again not accounted for because the molecular simulations show that only one layer can be formed [45] in such confined space. That will leave the Ising contribution W intact and creates a composite adsorption model. Despite the fact that Rutherford’s modification includes both sorbent–sorbate and sorbate–sorbate bonds, the one-layer limitations are obvious.

2.3. Compound Model

Looking at the CMMS model evolution and its solution at specific cases of its rate constants, it is apparent that the sigmoidal shape is difficult to achieve with a single model approach with just the classical approach of equivalent adsorption sites in the Langmuir model. The multisite adsorption concept is also used in the work of Do and Do [46,47] as a part of their combined adsorption model. Do and Do [46] assumed a more complex system and attempted to address the adsorption process from its individual parts; direct surface adsorption, adsorption in micropore space and capillary condensation in mesopore space. The original formula provided in the work of Do and Do [46] takes a form:

$$C_{\mu }=C_{\mu s}\frac{K_{\mu }{\sum }_{n=1}^i{x}^n}{K_{\mu }{\sum }_{n=1}^i{x}^n+K_{\mu }{\sum }_{n=1}^i{x}^{n-i}}+S_0\frac{K_f{\sum }_{n=1}n{x}^n}{1+K_f{\sum }_{n=1}{x}^n}$$

(15)

where $C_{\mu s}$ is the saturation concentration of the micropore, $K_{\mu }$ is the micropore equilibrium constant, $S_0$ is the functional group concentration and $K_f$ is the sorption equilibrium constant (chemisorption in the work of [46]). The first part of the equation describes microporous adsorption and the layer/cluster parameter was set to $n=6$ in the original work [46]. The meaning was to accommodate for forming of H₂O pentamers inside the microporous space as the minimum stable amount of molecules inside the hydrophobic environment as previously investigated by Kaneko et al. [48]. In their later work [47], the parameter n becomes adjustable and can be changed according to the fitting requirements. The second part of Equation (15) is then the well-known BET type equation describing the surface adsorption. Do and Do have realized that the problematic of capillary condensation in the mesoporous space cannot be easily described by the standard adsorption models and they suggested an alternative approach based on Kelvin equation for evaporation and Cohan equation for condensation. The capillary condensation section in the isotherm model is typically indicated by invariable energy on water uptake given the adsorption energy being proportional to the effective pore radius. This proportionality is given due to the change in surface energy caused by the curvature of the free liquid surface influencing the vapour pressure over the liquid surface as described by the Kelvin equation.

$$ln\frac{P}{P_0}=-\frac{2H\gamma V_l}{RT}$$

(16)

where $\frac{P}{P_0}$ is the relative pressure, H is the curvature of the meniscus in the capillary, $\gamma$ is the surface tension of the liquid nitrogen, the $V_l$ is the molar volume of liquid nitrogen. Therefore, the Kelvin equation can describe the dependency of critical pore radius $r_K$ (Kelvin radius) on the surface tension and contact angle based on its reciprocal relation to the meniscus curvature $r_K=1/H$. As a result that at the time of the condensation occurring a certain amount of the vapour is already adsorbed on the adsorbent’s walls, the final pore/capillary radius $r_P$ is calculated as a sum of the layer thickness $t_L$ and the critical radius $r_K$.

$$r_p=r_K+t_L$$

(17)

The layer thickness of already adsorbed vapour $t_L$ can be obtained by a statistical approach for adsorbates other than nitrogen. In the Do and Do work, the thickness is given for the activated carbon based on the standard adsorption model:

$$t_L=\frac{v_M{S}_0}{S_M}\frac{K_f{\sum }_{n}n{x}^n}{1+K_f{\sum }_n{n}^n}$$

(18)

where $v_M$ is the molar volume of the adsorbate in liquid form (condensed), and $S_M$ is the surface area of the mesopore space, while $S_0$ is again the concentration of the functional group.

2.4. Henry–Sips Combination Model

Another obvious, yet not often used, approach to the modelling of complex adsorption systems is a combination of the fundamental models to create hybrid models. This way, the differentiated parts of the adsorption process on a given adsorption pair can be described reasonably and more clearly. A combined model of Henry and Sips isotherms was proposed by Kim et al. [24]. The empirical combination of the original equations of Henry and Sips with a proportionality factor ${\beta }_{HS}$ is written in the following way:

$$q_e={\beta }_{HS}{K}_H\left(\frac{P}{P_0}\right)+\left(1-{\beta }_{HS}\right)\frac{q_0{\left(K_S\frac{P}{P_0}\right)}^{1/r}}{1+{\left(K_S\frac{P}{P_0}\right)}^{1/r}}$$

(19)

The right part of Equation (19) describes formally Langmuir–Freundlich isotherm model used by Sips [49], and it can be found in literature in several modifications [50]. In the above equation, the $q_e$ is the adsorption amount at the equilibrium, $q_0$ is the maximum adsorption amount (saturation uptake) and the $K_H$ and $K_S$ are Henry’s and Sips’ constants, respectively. The exponential factor r is then from the Langmuir–Freundlich model a measure of heterogeneity of the surface. The $\beta$ parameter is dependent on the partial pressure through an exponential giving a final constant ${\beta }_S$:

$${\beta }_{HS}=e^{-{\beta }_S\frac{P}{P_0}}$$

(20)

In this form, the beginning of the adsorption is governed by a site-controlled mechanism with limited adsorption at lower relative pressure. This combination was selected to describe two different processes of water adsorption on the example of ferroaluminophosphate zeolite AQSOA™ FAM-Z02. However, the fitting of this hybrid model is rather complicated and the absence of a strict boundary between the individual models makes the resulting model constants inconvenienced by superfluous errors.

2.5. Associated Energy Distribution Functions Model

An attempt to develop a universal isotherm model was further made by Ng et al. [18], building on a multi-site adsorption theory and site-dependent energy distribution functions. In this model, the heterogeneous surface is approximated according to the homotattic patch approximation (HPA). The HPA considers the heterogeneous surface as a formation of many homogenous patches, each patch described by its localised constants and connected by energy distribution function (EDF). The total adsorption uptake $\theta$ is then a sum of localised adsorptions ${\theta }_i$ as a function of energy $ϵ$ written as

$$\theta =\frac{N_a}{S_o}=\sum _{i=1}^{\infty }{\theta }_i\left(ϵ\right)$$

(21)

where $N_a$ and $S_o$ are then adsorbed molecules count and available adsorption sites number, respectively. Following the mathematical expansion of the sum of the local adsorption coverage based on their fractional availability, the sum of the energy distribution functions $X\left(ϵ\right)$ is proportional and equal to unity at the saturation point (endpoint of the adsorption):

$${\int }_0^{\infty }X\left(ϵ\right)dϵ=\sum _{n=1}^{\infty }\frac{s_i}{S_o}\left(ϵ\right)=1$$

(22)

where $s_i$ is the surface fractional availability (part of the surface) and $S_o$ is the overall adsorption area. In Equation (22) the total energy distribution function $X\left(ϵ\right)$ is formally a sum of respective energy functions corresponding to and defined by each localised adsorption site. Finally, the total adsorption uptake $\theta$ can be written as an integral of the total energy distribution over the whole system:

$$\theta ={\int }_0^{\infty }\left\{\theta \left(ϵ\right)X\left(ϵ\right)\right\}d\left(ϵ\right)$$

(23)

The adsorption model itself is then based on the Absolute Rate Theory [51], provided the localised adsorption uptake $\theta \left(ϵ\right)$ at any given time is equal to the difference in adsorption and desorption rate $R_{ad}$ and $R_d$, respectively:

$$\frac{d\theta \left(ϵ\right)}{dt}=R_{ad}-R_d$$

(24)

Therefore, according to the Absolute Rate Constant, the rate constants can be written based on the Arrhenius law as [52]:

$$\frac{d\theta \left(ϵ\right)}{dt}=K_{ad}P(1-\theta )exp\left(\frac{-{ϵ}_{ad}}{RT}\right)-K_d\theta exp\left(\frac{-{ϵ}_d}{RT}\right)$$

(25)

where P is the pressure, $K_{ad}$ and $K_d$ are respective pre-exponential constants and ${ϵ}_{ad}$ and ${ϵ}_d$ are the activation energies of adsorption and desorption, respectively. Naturally, at the equilibrium the uptake remains constant ($d\theta \left(ϵ\right)/dt=0$) and provided $K=K_{ad}/K_b$ and $\Delta ϵ={ϵ}_d-{ϵ}_{ad}$ we will get practically the Langmuir-type isotherm:

$$\theta \left(ϵ\right)=\frac{KPexp\left(\frac{\Delta ϵ}{RT}\right)}{1+KPexp\left(\frac{\Delta ϵ}{RT}\right)}$$

(26)

That means that only pairs site-particles are considered and the adsorption is finite. Furthermore, the model assumes that for accurate prediction of adsorption uptake is important the adsorption sites distribution over the surface of several heterogeneous patches. This approach is used frequently and the energy distribution approximations are used across the adsorption models [53]. Ng et al. [18] in their work suggest using symmetrical Gaussian function including fractional probability factor ${\gamma }_i$. Furthermore, because the probability factor denotes the share of the individual energy nodes on the whole system, the sum of all the factors has to be equal to unity $\sum {\gamma }_i=1$. Finally, the total energy distribution as the sum of the adsorption sites is defined by the individual Gaussian energy distributions as follows:

$$X\left(ϵ\right)=\sum _{i=1}^n{\gamma }_i{\left[\frac{exp\left(\frac{\Delta ϵ-{ϵ}_{oi}}{m_i}\right)}{m_i{\left(1+exp\left(\frac{\Delta ϵ-{ϵ}_{oi}}{m_i}\right)\right)}^2}\right]}_i$$

(27)

where ${ϵ}_{oi}$ is the adsorption energy with highest occurrence frequency which is in a practical sense the mean value shown as the peak of the energy distribution, and $m_i$ represents the energy variations depending on the heterogeneity of the adsorption site patch. Based on Equation (23), the final integrated model is then given as:

$$\theta =\sum _{i=1}^n{\gamma }_i{\left[\frac{{\left(\frac{P}{P_0}exp\left(\frac{{ϵ}_{oi}}{RT}\right)\right)}^{\frac{RT}{m_i}}}{1+{\left(\frac{P}{P_0}exp\left(\frac{{ϵ}_{oi}}{RT}\right)\right)}^{\frac{RT}{m_i}}}\right]}_i$$

(28)

Ng et al. [18] state that for isotherm Types I through V, the number of nominal peaks is two, n = 2, and for more complex systems, e.g., Type VI, it can be higher. This approach is, therefore, similar to the model based on statistical physics approach introduced by Khalfaoui et al. [54] and further modified by Yahia et al. [16]. The Ng’s model capabilities are also demonstrated on all isotherm types.

2.6. Specifics of Adsorption Process of Mesoporous Silica

Based on the general approach for approximation of the adsorption process described above, it seems that a statistical approach based on various modifications of adsorption distribution, or energy distribution, works very well for many systems with well-defined properties. However, many models of experimental data are often limited by only a partial fit due to many concurrent processes during adsorption. From the previous explanation, we can consider the adsorption process to consist of three main, energetically distinguishable, processes; monolayer adsorption, multilayer adsorption and capillary condensation.

In this work, only adsorption isotherms Type IV and Type V are considered based on the data of selected mesoporous silica samples TMPSs (Taiyo Kagaku Mesoporous Silica), SBA-15 and MCM-41 with water as the adsorbate. The shape and the complication of this particular adsorption pair is given by the fact that the adsorption process of the mesoporous materials is exactly the combination of two to three distinguishable processes depending on the surface characteristics as stated above. Freshly calcined mesoporous silica lacks silanol groups [55,56], because of the high preparation temperatures. The lack of nucleation sites and apolar surface consisting only from pure silica explains the Type V shape. However, after decent interaction and often elevated temperature with the water molecules, the silica surface is gradually hydrated, and silanol groups are formed on the surface [55,57]. This change creates a partially polar surface environment and enhances the multilayered adsorption of water vapour. The situation is different when a foreign atom with different electronegativity potential is introduced into the structure. Here, materials with the aluminium doped structure are used in comparison to aluminium-free materials to obtain both Type IV and Type V isotherm curves. The additional not atypical problem of isotherm modelling of these systems arises from the pore size, which can be occasionally on the boundary of the microporous and mesoporous region around 2 nm, especially when considering adsorption materials for energy conversion and conservation. That gives curves with a problematic region of influence between capillary condensation and multilayer adsorption in confined space as described by the CMMS or Do and Do models.

From the regular analysis and common approach to modelling, we can simulate any isotherm shape given enough independent constant in the model function. Comparable approach can be observed in models of Yahia et al. [16], Ng et al. [18] or Do and Do [46]. Although the approach is not flawed, the simulation and fitting of the experimental data can lead to an overly complex solution and without a necessary robust modelling method the simulations can give results which would be difficult to repeat or with limited physical meaning. Therefore, this work focuses on developing a reasonably simple approach that can provide highly accurate yet reproducible and repeatable results of adsorption model fitting for mesoporous materials especially those with developed capillary condensation region.

3. Materials and Methods

The experimental data in this work were obtained from industrially manufactured mesoporous silica materials from Taiyo Kagaku Co., Ltd. (Yokkaichi, Mie, Japan) and laboratory prepared samples of MCM-41 and SBA-15 mesoporous silica. The manufactured materials belong to the TMPS group and the selected five samples are specified by different pore size ranging from 1.8 $\mathrm{n}$$\mathrm{m}$ to 4.2 $\mathrm{n}$$\mathrm{m}$ and with or without aluminium doping. The aluminium doped samples are signed with a capital latter A at the end of the sample name. For a comparison, laboratory prepared MCM-41 and SBA-15 mesoporous silica were examined for the water adsorption as well. The MCM-41 was prepared by cationic surfactant CTMABr in ammonia solution as described by Grün et al. [58] and by adjusting the method to create a smaller pore size comparable to the studied TMPS specimens by decreasing the ammonia volume to 3 $\mathrm{m}$$\mathrm{L}$. The results are spherical cells with diameter spanning from 100 $\mathrm{n}$$\mathrm{m}$ up to 2 $\mathsf{\mu }$$\mathrm{m}$ on average with pore diameter approximately $3.3$ $\mathrm{n}$$\mathrm{m}$. The SBA-15 was synthesised according to Sayari et al. [59] with the Pluronic P123 template. To adjust the pore size, we chose the first micelle stabilisation phase to 20 h at $T=$ 35 ${}^{\circ }\mathrm{C}$ and the subsequent growing phase to $T=$ 100 ${}^{\circ }\mathrm{C}$ for 24 h. The result was standard rod-like cells with a diameter of 200$\mathrm{n}$$\mathrm{m}$–300$\mathrm{n}$$\mathrm{m}$ and length around 1 $\mathsf{\mu }$$\mathrm{m}$. The pore size was found out as $7.5$ $\mathrm{n}$$\mathrm{m}$. Similarly to the TMPS materials, both the SBA-15 and the MCM-41 were used without any modifications. The properties of all used materials are summarised in

Table 1. Materials properties; total surface area $S_m$ (BET), total pore volume $V_T$, mesopores pore volume $V_{mes}$, pore diameter $d_{DFT}$ and pore pitch $a_{XRD}$.

Sample	Total Surface Area ${\mathit{S}}_{\mathit{m}}$ ($\mathbf{m}$${}^2$ $\mathbf{g}$−1)	Total Pore Volume ${\mathit{V}}_{\mathit{T}}$ ($\mathbf{c}$$\mathbf{m}$${}^3$ $\mathbf{g}$−1)	Mesopores Pore Volume ${\mathit{V}}_{\mathbf{mes}}$ ($\mathbf{c}$$\mathbf{m}$${}^3$ $\mathbf{g}$−1)	Pore Diameter ${\mathit{d}}_{\mathbf{DFT}}$ ($\mathbf{n}$$\mathbf{m}$)	Pore Pitch ${\mathit{a}}_{\mathbf{XRD}}$ ($\mathbf{n}$$\mathbf{m}$)
TMPS-1.5	669.4	0.423	0.370	1.9	3.6
TMPS-1.5A	602.0	0.380	0.334	1.8	3.6
TMPS-2A	519.9	0.505	0.418	2.3	4.4
TMPS-4A	670.4	0.859	0.760	4.0	5.6
TMPS-4R	725.3	0.977	0.888	4.2	5.6
MCM-41	759.0	0.756	0.675	3.3	4.7
SBA-15	583.9	0.937	0.852	7.5	11.5

. The detailed investigation on the physical properties of the used mesoporous silica materials is provided elsewhere Mikšík et al. [60].

The adsorption isotherms were measured on high accuracy surface characterisation analyser 3Flex (Micromeritics, Norcross, GA, USA). The schematic of the measurement device is provided in Figure 2. The 3Flex analyser is a constant volume analyser and uses volumetric method to establish the adsorption amount. The standard measurement process follows a set routine based on sampling certain amount of the adsorbate vapour into the closed manifold space (closed by V1-6 valves) with known volume. In the second step the valve to the desired sample tube is opened (one of V1-3 valves) and the adsorbate is introduced to the sample. The sampling time is typically couple of seconds after which the valves is closed and the sampled vapour amount is calculated on the bases of the pressure difference before and after the introduction. The final adsorbed amount is then calculated based on introduced amount of adsorbed vapour and the residual pressure inside the sample tube after equilibrium is reached. The equilibrium time was set to 10 s with less than 0.01% change of the adsorption amount. The relative error of the measurement based on the reference material provided by the manufacturer was less than u < $0.01$. However, repeated measurements showed relative error of up to u = $0.05$ on the adsopriton amount and u = $0.02$ on the relative pressure as a joint error of weighting and vapour source preparation. However, with careful weight measurement the repeatability was typically 98. The relative deviation of the adsorption amount was also confirmed by an absolute thermogravimetric (TG) instrument MSB-TG (MicrotracBEL, Osaka, Japan) equipped with a magnetic suspension balance (MSB) (Rubotherm GmbH, Bochum, Germany). The high precision pressure transducer’s relative error was $u<0.0015$ of the reading accuracy allowing a negligible increase of the error through the measurement. The free volume of the experimental vessel and the sample was established through introducing a known amount of helium from the manifold to the fully evacuated sample tube. The volume of the sample tube containing the sample was then automatically calculated based on the pressure difference of the manifold and compared to the pressure increase inside of the sample tube before and after introducing the gas. Considering the size of the water molecules and the helium molecules, the deviation caused by the inaccessible space (very small crevices) is further on assumed as negligible.

Prior to the measurement, the samples were thoroughly degassed at 200 ${}^{\circ }\mathrm{C}$ and high vacuum ( $<0.01$ $\mathrm{Pa}$) for 5 h before weighting. Nitrogen gas was used to fill the manifold and sample tube to avoid unnecessary contamination during the manipulation between the measurements and weighting. After weighting, the samples were degassed once more at 200 ${}^{\circ }\mathrm{C}$ for another 5 h. The residual pressure over the samples in the sample cell was confirmed to be $<0.1$ $\mathrm{Pa}$ in all cases at the room temperature of 25 ${}^{\circ }\mathrm{C}$ at the end of the degassing process. The drying temperature of 200 ${}^{\circ }\mathrm{C}$ was selected with regards to the water adsorption properties of the possible presence of silanol (-Si-OH) groups [56]. Furthermore, to sustain high reproducibility, each measurement was conducted with a fresh sample every time. The TMPS samples were measured at five different temperatures $T=$ 5, 15, 25, 35 and 40${}^{\circ }\mathrm{C}$. The MCM-41 and SBA-15 were measured at four different temperatures 10, 20, 30 and 40${}^{\circ }\mathrm{C}$. The vapour source was prepared from ultra-pure deionised water with standard resistivity of $>18.2$ $\mathrm{M}$$\mathsf{\Omega }$$\mathrm{c}$$\mathrm{m}$. The vapour source was degassed by regular procedure employing liquid nitrogen and high vacuum source. The degassing cycle was repeated until the pressure of vapour source at the room temperature (25 ${}^{\circ }\mathrm{C}$) was equal to the referenced values (REFPROP).

3.1. New Type IV and Type V Adsorption Isotherm Model

Through the theoretical approaches to accommodate different adsorption processes described in the introduction in one model is deemed practically impossible. Especially challenging is the adsorption on mesoporous materials and those materials on the microporous and mesoporous boundary, as mentioned above as well. This problematic stands on the perspective of two synergic effects with different physical explanation and therefore governed by different properties. While the direct adsorption on the surface can be properly described by the standard adsorption theories such as Langmuir, BET or GAB, the limits of the capillary condensation in the mesoporous space lies in the mesopore volume distribution which roughly corresponds to the pore size distribution. As the condensation cannot be described simply by the empty adsorption sites nor through multilayer adsorption, as the limits are more of a geometrical constitution rather than chemical, a different approach should be considered for the mesoporous materials. In this work, an approach for the modelling was, therefore, taken from the perspective of the dependency of capillary condensation relative pressure directly on the pore size distribution. Such system was also briefly mentioned in the work of Do and Do [46] without further explanation per se. This approach, however, should require the identification of the pore size distribution beforehand or at least the adequate distribution model. Previously [60], the TMPS materials were analysed on their volumetric-based pore size distribution according to the Density Functional Theory (DFT) method, and a surprisingly good fit was observed with a modified standard distribution function:

$$y\left(x\right)=\alpha \frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{\frac{-{\left(x-\omega \right)}^2}{2{\sigma }^2}}$$

(29)

where x is the pore size, $\alpha$ is a dimensional coefficient and $\sigma$ and $\omega$ are regular coefficients for the standard distribution function, where $\sigma$ correlates to the range and $\omega$ is the position factor. Considering that the adsorption process for mesoporous materials is a combination of adsorption and capillary condensation, we can write:

$$q_e=q_{ea}+q_{ec}$$

(30)

where the adsorption amount at equilibrium $q_e$ is a sum of the regular adsorption $q_{ea}$ and capillary condensation $q_{ec}$. By further developing this dependency it is possible to deduce a combination of standard adsorption model and condensation process based on the pore distribution. Firstly, it is necessary to rewrite the Equation (29) into an integral form while $y\left(x\right)=q_{ec}\left(\frac{P}{P_0}\right)$:

$$\frac{d{q}_{ec}}{d\frac{P}{P_0}}=\alpha \frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{\frac{-{\left(\frac{P}{P_0}-\omega \right)}^2}{2{\sigma }^2}}$$

(31)

The boundary conditions of the Equation (31) are logically stated for the states in which the isotherm is defined from $P/P_0$ = 0 to $P/P_0$ = $P/P_0$, while the upper limit is unity at full system saturation where $P=P_0$. The integration of this modified standard distribution then leads naturally to error function and after some algebraic manipulation, we will get the cumulative condensation amount:

$$q_{ec}=\alpha \frac{1}{2}\left[1+erf\left(\frac{\frac{P}{P_0}-\omega }{\sqrt{2}\sigma }\right)\right]$$

(32)

As the adsorption on hydrophilic materials is multilayer in principle, the selection of the surface adsorption part was obvious with respect to the BET theory. We selected two well-known models based on the assumption of the practicality of the new model; the multilayer BET (Equation (7)) and the modified BET model known as the GAB model. The main difference between these models would be, therefore, the expectation of the finite number of adsorbed layers (i.e., confined space) in opposition to the reduced relative pressure $x=K_a\frac{P}{P_0}$ to compensate for the saturation pressure deviation in case of the GAB model. The multilayer BET model extended by the mesopore condensation distribution function defined by Equation (32) and the final equation takes the following form (further on only as mBET+G):

$$\begin{array}{cc}\hfill q_e& =q_m\frac{Cx}{1-x}\frac{1-\left(n+1\right)x^n+n{x}^{n+1}}{1+\left(C-1\right)x-C{x}^{n+1}}\hfill \\ \hfill & +\alpha \frac{1}{2}\left[1+erf\left(\frac{\frac{P}{P_0}-\omega }{\sqrt{2}\sigma }\right)\right]\hfill \end{array}$$

(33)

where the variable x is the standard relative pressure $x=P/P_0$. Above that, we can describe the individual constants from Equation (33) by the temperature dimension T, which will be demonstrated in the following parts.

The combination of the condensation distribution function of Equation (32) part with the GAB model will take a similar route giving the final form:

$$\begin{array}{cc}\hfill q_e& =q_m\frac{C_{G}x}{\left(1-x\right)\left(1+\left(C_G-1\right)x\right)}\hfill \\ \hfill & +\alpha \frac{1}{2}\left[1+erf\left(\frac{\frac{P}{P_0}-\omega }{\sqrt{2}\sigma }\right)\right]\hfill \end{array}$$

(34)

In Equation (34) (further on referred as GAB+G) x is the standard relative pressure substituted by the reduced relative pressure from Equation (8). As mentioned before, the two new models have two parts describing individual processes during the adsorption. The GAB modification has been proven as a valid model for many systems surpassing the original BET model especially at higher relative pressures; however, it is derived from the assumption of infinite adsorption layers, which can lead to significant overestimation in adsorption. On the other hand, it is simpler than the second model based on the finite-layer derivation of the BET theory. Similarly to BET model the $C_G$ constant represents the affinity of the adsorbate to the adsorbent.

The $\alpha$ and $q_m$ (first layer adsorbed volume), are properties of the system and should be independent on the temperature. However, in the case of $q_m$, it is more appropriate to talk about apparent first layer capacity, as it is not necessarily an exact value, but rather an approximation of the model itself.

3.2. Modelling, Error Analysis and Goodness of Fit

To assess the quality of the model and the goodness of fit, several methods have been employed. As the standard method for the goodness of fit was used the non-linear Chi-Square Test (${\chi }^2$) in the following form:

$$\left({\chi }^2\right)=\sum _{i=1}^n\frac{{\left(q_{cl}-q_{ex}\right)}^2}{q_{ex}}$$

(35)

where $q_{cl}$ is the calculated value and $q_{ex}$ is the experimentally measured value. (${\chi }^2$) method also has its shortcomings as it is supposed to be fit to data with a certain distribution. In case that the model fits well and the distribution of the points is too small, the results are small as well and a direct comparison can prove to be impractical. Hence, for standard comparison of dependencies, the usual coefficient of determination $R^2$ is used as well to assess the correlation between the models and experimental data. The equation is used in the form:

$$R^2=1-\frac{\sum {\left(q_{cl}-q_{ex}\right)}^2}{\sum {\left(q_{cl}-{\overline{q}}_{ex}\right)}^2}$$

(36)

where the parameters are the same as in Equation (35), while ${\overline{q}}_{ex}$ is the average value of the adsorbed amount. As a result that the coefficient of determination does not provide adequate information when the data range spans through several orders of magnitude, to assess the absolute fit of the calculated model on the experimental data an average absolute deviation (AAD) and average relative deviation (ARD%) were used as a direct measure for the data fitting. The ARD% was used in its natural form:

$$ARD\%=\frac{{\displaystyle \sum _{i=1}^n\frac{{\displaystyle ∥q_{cl}-q_{ex}∥}}{{\displaystyle q_{ex}}}}}{{\displaystyle n}}$$

(37)

The average absolute deviation was then calculated similarly without the scaling normalization by the effective experimental data $q_{ex}$:

$$AAD=\frac{{\displaystyle \sum _{i=1}^n∥q_{cl}-q_{ex}∥}}{{\displaystyle n}}$$

(38)

The parameters of the function are again identical to Equation (35) of the ${\chi }^2$ test.

The model simulation was conducted through a repeated non-linear regression fitting using Trust Region Reflective (TRR) algorithm. An important criterion was an introduction of a weighting factor to compensate for the scalability of the adsorption amount at low relative pressure with optimization of the absolute average deviation fitting parameter. The fundamental fitting (preliminary fitting) of the models was executed with a different constant weighting factor according to $ln\left(\frac{P}{P_0}\right)/N$ distribution, where N is the number of the data points. This method was selected to better compensate for the imbalance in the input data, giving the beginning of the adsorption higher priority on a more local level. This way, a good fit for the basic model could be ascertained, providing us with meaningful results in terms of the BET adsorption theory and very good fit with experimental data. The same approach was applied to the other models as well. All the calculations and modelling of the isotherm curves were executed by original Python scripts.

4. Results

From the explanation above, it is evident that simple models would be difficult to apply to mesoporous materials as materials with strong duality (or triality) in their adsorption mechanisms. For the comparison of the new isotherm model, we have chosen the Rutherford model [42], Henry–Sips model [24] and Ng’s model [18] with two nominal peaks ($n=2$) (Ng-2f) as suggested in their work for the isotherm types IV and V. The selection of the comparison models was done based on the best fit out of the regular models described previously in this work. Although the goodness of fit of these models somehow lacks certain aspects, the results of the selected models prove to provide sufficient fit with the present mesoporous silica materials and the results are comparable to the new models. A full summary of the final form, at which all the models were used, is provided in

Table 2. Summary of the models in their final form for the comparative simulations.

Model Name	Label	Equation	Reference
Henry–Sips model	Henry–Sips	$q_e={\beta }_{HS}{K}_H\left(\frac{P}{P_0}\right)+\left(1-{\beta }_{HS}\right)\frac{q_0{\left(K_S\frac{P}{P_0}\right)}^{1/r}}{1+{\left(K_S\frac{P}{P_0}\right)}^{1/r}}$	Ref. [24]
Ng’s two part model	Ng-2f	$q_e=q_0{\gamma }_1\frac{{\left(\frac{P}{P_0}exp\left(\frac{{ϵ}_{o1}}{RT}\right)\right)}^{\frac{RT}{m_1}}}{1+{\left(\frac{P}{P_0}exp\left(\frac{{ϵ}_{o1}}{RT}\right)\right)}^{\frac{RT}{m_1}}}+q_0{\gamma }_2\frac{{\left(\frac{P}{P_0}exp\left(\frac{{ϵ}_{o2}}{RT}\right)\right)}^{\frac{RT}{m_2}}}{1+{\left(\frac{P}{P_0}exp\left(\frac{{ϵ}_{o2}}{RT}\right)\right)}^{\frac{RT}{m_2}}}$	Ref. [18]
Rutherford’s model	Rutherford	$q_e=\frac{q_{satL}{b}_L\frac{P}{P_0}}{1+b_L\frac{P}{P_0}}+\frac{q_{sat}{K}_0\frac{P}{P_0}}{K_0\frac{P}{P_0}+W^2}$ , $W=\frac{1}{2}\left[1-K_1\frac{P}{P_0}+\sqrt{{\left(1-K_1\frac{P}{P_0}\right)}^2+4K_0\frac{P}{P_0}}\right]$	Ref. [43]
Multilayer BET + modified distribution function model	mBET+G	$q_e=q_m\frac{Cx}{1-x}\frac{1-\left(n+1\right)x^n+n{x}^{n+1}}{1+\left(C-1\right)x-C{x}^{n+1}}+\alpha \frac{1}{2}\left[1+erf\left(\frac{\frac{P}{P_0}-\omega }{\sqrt{2}\sigma }\right)\right]$, $x=\frac{P}{P_0}$	This work
GAB + modified distribution function model	GAB+G	$q_e=q_m\frac{C_{G}x}{\left(1-x\right)\left(1+\left(C_G-1\right)x\right)}+\alpha \frac{1}{2}\left[1+erf\left(\frac{\frac{P}{P_0}-\omega }{\sqrt{2}\sigma }\right)\right]$, $x=K_a\frac{P}{P_0}$	This work

.

The graphical representation of the new mBET+G and GAB+G models compared to the experimental data and literature models is provided in Figure 3 on the example of TMPS-2A and SBA-15. It is visible that all the selected models perform well except for the Henry–Sips model which cannot properly describe and follow the two modes of the adsorption process which is even more apparent given the respective examples in Figure 3 which are representatives of both Type V (SBA-15) and Type IV (TMPS-2A) isotherms. To objectively compare the fitting of the models, a numerical comparison is provided by the means of ${\chi }^2$ and ARD% for all tested samples. The statistical parameters are summarized in

Table 3. Statistical coefficients overview (ARD% and ${\chi }^2$) of Henry–Sips, Ng-2f, Rutherford, mBET+G and GAB+G models for all samples at 25 ${}^{\circ }\mathrm{C}$ (20 ${}^{\circ }\mathrm{C}$ for SBA-15 and MCM-41).

Sample	Henry-Sips		Ng-2f		Rutherford		mBET+G		GAB+G
	${\mathit{\chi }}^2$	ARD%	${\mathit{\chi }}^2$	ARD%	${\mathit{\chi }}^2$	ARD%	${\mathit{\chi }}^2$	ARD%	${\mathit{\chi }}^2$	ARD%
TMPS-1.5	0.00966	9.4	0.00516	7.4	0.00557	8.1	0.00484	7.7	0.00522	8.0
TMPS-1.5A	0.05253	20.2	0.01614	21.1	0.03204	16.0	0.02959	15.2	0.03135	15.6
TMPS-2A	0.09132	19.1	0.03296	16.5	0.05810	15.8	0.05327	15.3	0.05942	16.0
TMPS-4A	0.08526	18.9	0.04377	12.5	0.05490	15.3	0.03431	12.5	0.04811	14.7
TMPS-4R	0.05619	10.9	0.04911	13.9	0.02227	8.5	0.03120	10.5	0.03413	11.0
SBA-15	0.81700	16.6	0.31625	15.4	0.24631	9.3	0.04780	4.1	0.17025	11.1
MCM-41	0.04402	11.2	0.01194	10.6	0.00868	8.7	0.00975	9.3	0.01261	10.9

.

The best fit from the regular models based on the ${\chi }^2$ and ARD% values shows Rutherford’s model following the CMMS approach; moreover, this model works with only five constants making it more statistically significant. The visible improvement of the fitting is apparent from the ARD% values, which had been decreased in almost all cases. The better fitting can also be observed directly from the graphical output of the residual analysis in Figure 4. However, when comparing the individual models, we have to evaluate all of the statistical factors.

Even if the relative deviations are rather high (mostly > 10%) the average absolute deviations listed in the

Table 4. Adsorption models’ fitting parameters and constants for all samples at 25 ${}^{\circ }\mathrm{C}$ (20 ${}^{\circ }\mathrm{C}$ for SBA-15 and MCM-41) with average absolute deviation (AAD) of the fit of the experimental data.

Henry-Sips
Sample	${\mathit{\beta }}_{\mathit{S}}$	${\mathit{K}}_{\mathit{H}}$	${\mathit{K}}_{\mathit{S}}$	$\mathit{r}$	${\mathit{q}}_{\mathbf{0}}$		AAD
TMPS-1.5	2.727	0.7366	3.305	0.06731	0.3077		0.00452
TMPS-1.5A	8.080	1.8246	3.635	0.16267	0.3156		0.00709
TMPS-2A	2.684	0.9713	2.586	0.07013	0.3721		0.00780
TMPS-4A	1.114	0.5472	1.594	0.02451	0.9419		0.01035
TMPS-4R	1.000	0.3322	1.457	0.00884	1.2734		0.01431
SBA-15	1.000	0.7240	1.237	0.00917	1.0224		0.02725
MCM-41	1.934	0.5992	1.934	0.02348	0.6673		0.01034
Ng-2f
Sample	$\mathit{\alpha }$	${\mathit{ϵ}}_{\mathbf{1}}$	${\mathit{m}}_{\mathbf{1}}$	${\mathit{ϵ}}_{\mathbf{2}}$	${\mathit{m}}_{\mathbf{2}}$	${\mathit{q}}_{\mathbf{0}}$	AAD
TMPS-1.5	0.563	2527.64	2636.23	2897.43	164.25	0.4182	0.00285
TMPS-1.5A	0.692	573.21	4228.63	3153.63	208.55	0.5029	0.00366
TMPS-2A	0.647	1.08	3839.92	2325.03	162.66	0.6511	0.00518
TMPS-4A	0.435	1.16	3220.25	1138.37	68.87	0.9485	0.00636
TMPS-4R	0.344	7.38	1554.36	926.28	22.72	1.046	0.01332
SBA-15	0.608	1.49	1342.89	512.11	17.71	1.2928	0.01627
MCM-41	0.518	1582.88	56.72	1.11	2429.77	0.85	0.00469
Rutherford
Sample	${\mathit{b}}_{\mathit{L}}$	${\mathit{K}}_{\mathbf{0}}$	${\mathit{K}}_{\mathbf{1}}$	${\mathit{q}}_{\mathit{satL}}$	${\mathit{q}}_{\mathit{sat}}$		AAD
TMPS-1.5	4.48	0.01520	3.202	0.1718	0.2061		0.00289
TMPS-1.5A	35.35	0.06150	3.571	0.0927	0.2350		0.00523
TMPS-2A	14.57	0.01840	2.553	0.1342	0.3028		0.00664
TMPS-4A	4.41	0.00120	1.583	0.2016	0.5778		0.00609
TMPS-4R	0.31	0.00010	1.453	0.6087	0.7270		0.01175
SBA-15	0.48	0.00009	1.234	1.0000	0.5696		0.01819
MCM-41	1.68	0.00100	1.914	0.2661	0.4707		0.00470
mBET+G
Sample	$\mathit{C}$	$\mathit{n}$	$\mathit{\alpha }$	$\mathit{\sigma }$	$\mathit{\omega }$	${\mathit{q}}_{\mathit{m}}$	AAD
TMPS-1.5	7.32	2.22	0.1854	0.03552	0.3102	0.1111	0.00246
TMPS-1.5A	46.51	2.641	0.1896	0.05181	0.2734	0.0802	0.00461
TMPS-2A	20.97	3.479	0.2238	0.04155	0.3877	0.1035	0.00553
TMPS-4A	19.59	7.044	0.4666	0.02492	0.6305	0.0831	0.00537
TMPS-4R	2.58	7.415	0.6547	0.00994	0.6879	0.06	0.01231
SBA-15	6.25	13.809	0.4006	0.00753	0.8106	0.0859	0.01096
MCM-41	4.01	3.369	0.4268	0.01949	0.5221	0.1095	0.00416
GAB+G
Sample	${\mathit{C}}_{\mathit{G}}$	${\mathit{K}}_{\mathit{a}}$	$\mathit{\alpha }$	$\mathit{\sigma }$	$\mathit{\omega }$	${\mathit{q}}_{\mathit{m}}$	AAD
TMPS-1.5	35.59	0.1294	0.1892	0.03619	0.3107	0.169	0.00246
TMPS-1.5A	97.93	0.3969	0.1997	0.05393	0.2729	0.0888	0.00452
TMPS-2A	32.21	0.448	0.2384	0.044	0.3875	0.1295	0.00559
TMPS-4A	14.66	0.6896	0.4805	0.02596	0.6296	0.1132	0.00723
TMPS-4R	0.2	0.2902	0.6537	0.00989	0.6879	2.3111	0.01234
SBA-15	0.2	0.5105	0.4009	0.00717	0.8108	1.7683	0.01389
MCM-41	11.15	0.1	0.4297	0.01964	0.5222	0.3591	0.00432

, together with the full summary of the models’ constants parameters, are relatively small for all the models. The implications are clear when we directly confront these two values. If the AAD is small, the overall fit is acceptable, however, if the ARD% value is also small, even the fitting at lower relative pressures fits the experimental data. Then again, with ARD% values too high, the fit at the lower relative pressures is typically poor as it is apparent in Figure 5 showing detail fitting in the relative pressure $P/P_0$ region of 0 to 0.2. This comparison is indeed crucial if fitting the whole adsorption pressure range by physically meaningful parameters of the used models.

The new models are easily comparable to the recent models and in several aspects, they outperform them. The direct comparison based on ARD% reveals that far best results are reached by the mBET+G, which reaches the best ARD% values and ${\chi }^2$ parameters on average. That is, some of the materials show better fitting parameters for different individual models, however, the overall universal applicability is in favour of the mBET+G model. As an example can serve the material TMPS-4R which fits the mBET+G model with ARD% = $10.5$%, however, Rutherford’s model shows ARD% = $8.5$%, which is clearly better value. Similar results are also in case of the absolute deviation AAD of the TMPS-4R where Rutherford’s model performs the best among all the models with AAD = $0.0118$ $\mathrm{g}$ $\mathrm{g}$⁻¹ while the mBET+G is the second best with AAD = $0.0123$ $\mathrm{g}$ $\mathrm{g}$⁻¹. Afterwards, the mBET+G model shows better results and fitting parameters in all other materials compared to Rutherford’s model. The same situation is in case of TMPS-1.5A where the best performing model is Ng’s two-part model (Ng-2f) reaching very low AAD = $0.0037$ $\mathrm{g}$ $\mathrm{g}$⁻¹. The mBET+G and GAB+G then show slightly higher values of AAD = $0.0046$ $\mathrm{g}$ $\mathrm{g}$⁻¹ and AAD = $0.0045$ $\mathrm{g}$ $\mathrm{g}$⁻¹, respectively. The fitting example of the TMPS-4R and TMPS-1.5A is available in Figure 6. The results are therefore strongly in favour of mBET+G and the GAB+G model is used mainly for demonstration and comparison.

A strong point of the proposed combination is the complementarity of the two models; multilayer BET model and capillary condensation distribution function model. Comparing the fitting parameters of the individual models, the newly proposed model systems, especially the mBET+G, provide the most regular parameters in relation to the properties of the adsorbents. Looking at the Henry–Sips model, the proportionality factor ${\beta }_{HS}$ almost completely separates the individual parts of Henry contribution and Sips contribution at a certain point. That is not necessarily wrong for many of the materials which exhibit a strong Henry region, however, the applicability is strongly limited and, as shown in this work, the fitting is mostly insufficient. Different situation is in case of the Ng-2f’s model where the energy distribution functions are designed to overlap each other to create the resulting isotherm. Even though that the fitting is improved, the fitting parameters are not following the individual adsorption phases. That is apparent from irregular values of the parameters (very high or very small values of ${ϵ}_1,m_1,{ϵ}_2,m_2$) for different adsorbents regardless the properties. The best fit in comparison to the proposed mBET+G, therefore, provides the Rutherford’s model. Rutherford’s model is technically a combination of Langmuir isotherm and modified (simplified) CMMS model. However, it is important to note that Rutherford’s model is not primarily intended to cover the capillary condensation as was explained in the Section 2.2. Thus, the parameters of the model have limited meaning for materials with capillary condensation. As a result of this, the Rutherford’s model tends to overestimate the adsorption quantity at lower relative pressure for materials with low affinity to the adsorbate (SBA-15) or underestimate it for materials with high affinity (TMPS-4A) as clearly visible in Figure 5. From this perspective is clear that the proposed model has superior properties when considering materials with certain amount of capillary condensation. It would be possible to argue that a very good fit can be reached by expanding the Ng’s multipart model [18] or Yahia’s multipart model [16], however, this fitting would be achieved on the expenses of increased amount of necessary parameters and complexity.

From the Figure 5 and Figure 6 on the examples TMPS-4A and TMPS-1.5A, it is notable that most of the models underestimate the adsorbed quantity at low $P/P_0$, and therefore, underestimates the adsorption potential. This situation is most apparent in the case of material with high affinity to water and strong adsorbent–adsorbate interaction of the adsorption isotherm Type IV. Furthermore, at the end of the adsorption at high relative pressure, all of the samples show a continuous increase in the adsorption amount associated with progressing multilayer adsorption and at the further end with inter-particle condensation. This part is also very well described by the extended BET combined models as the multilayer formation is usually considered to continue through the whole pressure span. However, because the BET theory does not account for the inter-particle condensation, it may lead to inevitable compensation error, as shown in detail in Figure 4. The same problems, naturally, concerns other models as well as they do not explicitly deal with the physical properties of the system such as the shape and size of the adsorbent particles. Nevertheless, this region is seldom a concern during modelling or simulation of the adsorption systems and it is primarily used for the material analysis only.

Temperature Dependent Fitting and Analysis

For the purpose of a combination with a simulated environment, the model parameters should be fitted into a temperature function series. Hence, the creation of a universal model applicable for the whole temperature region was attempted in this work as well. Considering the best fit of the experimental data and the purpose of this study, we have selected only the mBET+G model for this procedure as it gives better overall performance over the GAB+G model. The approximation of the model’s parameters dependence on the temperature can be expressed as a linear function of temperature:

$$f\left(T\right)=p_{A}T+p_B$$

(39)

where T is the variable temperature, and $p_A$ and $p_B$ are the individual parameter’s temperature change rate and temperature-independent constant component, respectively. The linear approximation was selected for its simplicity and with the focus on sustainable accuracy which was validated in comparison for the individual isothermal data fitting. However, we need to emphasise that the linear approximation used specifically in this temperature range will not be necessarily also applicable beyond the scope of the presented experimental data. Hence, extrapolation of the data outside of the experimentally acquired temperature range should be validated.

The full approximation can be simply made based on the results of the individual isotherm analyses and their linear estimation. That gives, however, rather crude results with fitting errors far surpassing the original fitting parameters and a further optimization has to be applied. With the fundamental estimation taken from the simple linear regression of the individual parameters, the whole model was optimized on the whole data set at once through a continuous non-linear regression using the ARD% results as the weighting factors. This approach allowed us to have good agreement with the original fitting focusing on the sustainable fit at lower relative pressure and at the same time, we could partially compensate for the differences of the separate isotherm measurements. The final and optimized fitting is presented in Figure 7 on the example of TMPS-1.5 and TMPS-2A, which have different surface characteristics and mean pore size. The average relative deviation is typically higher than in case of the individual isotherm modelling and reaches AAD$\%=19.1$% for TMPS-2A compared to AAD$\%=15.3$% of the one isotherm fitting at 25 ${}^{\circ }\mathrm{C}$. Somewhat atypical is the fitting of TMPS-1.5 which shows overall deviation decrease in the temperature approximation over the one isotherm fitting dropping from the original AAD$\%=7.7$% at 25 ${}^{\circ }\mathrm{C}$ to the overall AAD$\%=7.4$%. This is caused by the slightly deficient fit at 25 ${}^{\circ }\mathrm{C}$ while the average of all isotherms of TMPS-1.5 was in fact $6.9$%. The fitting parameters of the temperature approximation for the whole sample set are summarized in

Table 5. mBET+G model’s temperature dependent parameters approximation. Valid temperature region 5 ${}^{\circ }\mathrm{C}$ to 40 ${}^{\circ }\mathrm{C}$.

Sample	${\mathit{C}}_{\mathit{A}}$	${\mathit{C}}_{\mathit{B}}$	${\mathit{n}}_{\mathit{A}}$	${\mathit{n}}_{\mathit{B}}$	${\mathit{\alpha }}_{\mathit{A}}$	${\mathit{\alpha }}_{\mathit{B}}$	${\mathit{\sigma }}_{\mathit{A}}$	${\mathit{\sigma }}_{\mathit{B}}$	${\mathit{\omega }}_{\mathit{A}}$	${\mathit{\omega }}_{\mathit{B}}$	${\mathit{q}}_{\mathit{mA}}$	${\mathit{q}}_{\mathit{mB}}$	${\mathit{\chi }}^{\mathbf{2}}$	ARD%	AAD ($\mathbf{g} {\mathbf{g}}^{-1}$)
TMPS-1.5	−0.1473	10.59	−0.0121	2.423	−0.00031	0.1940	−0.000094	0.0379	0.001235	0.2845	0.000750	0.09741	0.0427	7.4	0.0037
TMPS-1.5A	−0.7607	63.08	−0.0051	2.773	−0.00056	0.2036	−0.000090	0.0544	0.001087	0.2534	0.000307	0.07319	0.2225	19.8	0.0051
TMPS-2A	−0.3433	28.34	−0.0077	3.695	−0.00049	0.2372	0.000054	0.0393	0.001163	0.3634	0.000201	0.09833	0.3404	19.1	0.0065
TMPS-4A	−0.3528	31.59	−0.0420	8.461	0.00014	0.4556	0.000004	0.0253	0.000652	0.6223	0.000073	0.07974	0.4397	16.7	0.0169
TMPS-4R	−0.1422	7.60	−0.1538	13.176	0.00069	0.6223	0.000011	0.0116	−0.000264	0.7069	0.000485	0.04090	3.8124	21.3	0.0587
SBA-15	−0.0731	6.92	0.1291	9.647	−0.00111	0.4163	−0.000179	0.0129	−0.000259	0.8143	−0.000061	0.09365	0.7643	6.7	0.0225
MCM-41	−0.0010	4.11	0.0999	0.678	−0.00251	0.5088	−0.000523	0.0427	−0.000354	0.5472	−0.000949	0.12628	2.5692	19.5	0.0427

.

Based on the linear expression, we can distinguish two separate patterns within the parameters. The first group describing the adsorbent–adsorbate interactions consists of highly and partially temperature-dependent parameters C and n. These parameters decrease with temperature in all cases, which follows the general adsorption theory that higher adsorption temperature negatively influences the adsorbent–adsorbate adsorption energy resulting in weaker bonding. The second group consists of parameters which are directly dependent on the physical properties of the adsorbent and mostly influences by the adsorbate–adsorbate interactions, such as capillary condensation. Into this group also belongs the maximum surface uptake $q_m$ as it does describe the unchanging surface area of the adsorbent even though it is affected by the adsorbate-adsorbent to a certain degree. The values of $q_{mA}$ for most of the tested materials are well below $0.001$ $\mathrm{g}$ $\mathrm{g}$⁻¹ ${}^{\circ }\mathrm{C}$⁻¹ of the slope making it practically constant in the selected temperature range. Similar observations were noted for the parameters describing the capillary condensation region in the form of the distribution function consisting of $\alpha$, $\sigma$ and $\omega$ parameters, here with respective linear coefficients denoted as sub ${}_A$ and intercept constant denoted as sub ${}_B$. If we look closely to the effective values of those parameters, we can see that the linear coefficients are typically two to three orders smaller than the intercept constants. This mean that within the selected temperature range 5 ${}^{\circ }\mathrm{C}$ to 40 ${}^{\circ }\mathrm{C}$ the temperature dependence is limited to less than 10% on the final values. Furthermore, individual changes can be explained by the governing processes during the adsorption giving them physical meaning, although the explanation is not that straightforward as we would expect. For example, the change with temperature of the $\omega$ parameter, which describes the mean value of the capillary condensation relative pressure, is connected to the change of the surface tension from the Kelvin Equation (16). The increase in the surface tension caused by the temperature rise also causes the increase of the Kelvin radius leading to shifting of the condensation in the pores to higher relative pressure values. However, at the same, the density of the adsorbed layers decreases, causing an opposing force to this process as the adsorbed layer thickness increases on the surface of the adsorbent. Although this effect may be small in comparison, many of these processes affect the final value requiring further insight.

5. Discussion

Despite the similar origins, the materials presented in this work have different properties and behaviour arising from their chemical and geometrical structure. The typical and pure mesoporous silica usually exhibit limited adsorption at lower relative pressures as shown in Figure 3 on the typical representative SBA-15. A similar situation is in case of TMPS-4R and MCM-41, which have unaltered simple silica structure. These low hydrophilic properties of their surface are adequately described by the low C constants values from the mBET+G and GAB+G models. The C constant in the BET theory describes the intensity of the adsorbent–adsorbate interaction, and it is typically given by the following formula:

$$C=exp\left(\frac{E_1-E_L}{RT}\right)$$

(40)

where $E_1$ and $E_L$ are the average adsorption heat of the first layer and heat of liquefaction, respectively. The gas constant R is then used together with the thermodynamic temperature T. However, even though that the C and $C_G$ constants of GAB and BET models share the same origins, the values significantly vary as shown in Table 4 when the GAB+G model is showing considerably higher values, even several times higher than the mBET+G model. If we look further into the meaning of the C constant and its relation to other constants in the GAB and BET models as was explained before, we come to the realisation that in case of the GAB model the $C_G$ constant is practically co-dependent on the $K_a$ constant while in the BET model the relation is clearer from the original expression of Equation (7) where $x=P/P_0$. That means while in the GAB+G model the $C_G$ constant requires a more careful approach, in case of the mBET+G model Equation (40) can be applied directly provided the capillary condensation expression part does not interfere with the BET contribution to the model. Truly, in Figure 8 is clearly visible that the distribution function applies only in a very narrow window corresponding neatly to the capillary condensation region. At the same time, the multilayer adsorption following the BET model is undisturbed and continues to add to the adsorbed quantity. A certain flaw in this model can be perceived the sudden change in the surface area after the full pore filling making the inner surface of the pores completely inaccessible which is not accounted for in either part of the equation. However, despite that, the model continues without significant deviation from the experimental data until high relative pressure where the inter-particle condensation takes place. The inter-particle condensation is characterised by a sudden rise in the adsorption quantity beyond the $P/P_0>0.9$ mark and it depends on the particle size and shape. As a result that the inter-particle condensation is of low importance in the practical meaning, we tried to avoid this region in the experimental setup actively and only trace evidence of this process can be observed in Figure 6 on the example of TMPS-4R.

In the case of the aluminium doped samples TMPS-1.5A, TMPS-2A and TMPS-4A, the C values are significantly higher, showing strong hydrophilic properties. As a result that the aluminium content is similar to all of the augmented samples, the C should be as well similar. Yet, in fact only the TMPS-2A and TMPS-4A show similar values of $C=20.97$ and $C=19.59$, respectively, and the TMPS-1.5A shows C value as high as $C=46.51$. This significant increase is caused by the boundary conditions of the microporous region to which the TMPS-1.5A belongs with the average pore diameter of $d_{DFT}=1.8\mathrm{n}\mathrm{m}$. The closeness of the microporous region also affects the slightly higher value of the TMPS-2A compared to the TMPS-4A. The microporous conditions then also explain the higher values of the TMPS-1.5 compared to the other aluminium-free samples. The seemingly higher value of SBA-15 at $C=6.25$ is similarly caused by a certain content of microporous space created by the inter-pore channels as the remnants of the templating [61] and as found out as well during the detailed analysis of the sample characteristics [60].

If we attempt to compare the calculated values of the adsorption heat from Equation (40) of the first layer with the calculated values from experimental data by Clausius–Clapeyron (CC) relation from our previous work [26] the average values appear to be in agreement with those calculations. The highest value of the adsorption heat for the first layer based on the mBET+G is reached by the TMPS-1.5A with $E_1=55.4\mathrm{k}\mathrm{J} \mathrm{mol}$⁻¹, while the lowest one is TMPS-4R with $E_1=48.2\mathrm{k}\mathrm{J} \mathrm{mol}$⁻¹. The equivalent values from the Clausius-Clapeyron relation [26] for the first layer (calculated as an average value from $q=0\mathrm{mol} \mathrm{g}$⁻¹ to $q=q_m$) are at $\Delta H_{ads}=49.1\mathrm{k}\mathrm{J} \mathrm{mol}$⁻¹ and $\Delta H_{ads}=44.7\mathrm{k}\mathrm{J} \mathrm{mol}$⁻¹ for TMPS-1.5A and TMPS-4R, respectively. Considering that the calculated values of the adsorption heat by the CC relation are typically lower than the actual experimental results [62], the difference between these two samples is similar and equivalent to their respective surface characteristics.

From the fitting parameters of the mBET+G model is evident that the normal distribution works very well for all of the studied material samples. However, the normal distribution function used to describe the capillary condensation is still somehow lacking in terms of the physical meaning of the individual constants. Therefore, in here, we attempt to explain the individual constant using the physicochemical parameters, whereas the comparison in this work was possible thanks to the variety of the samples, varying by the pore size and surface characteristics. While the physical meaning of the C and n constants from the original BET model is clear and was explained above, parameters describing the capillary condensation may have a broader meaning.

The function describing the capillary condensation is built around the physical pore size distribution, which proved to the adjusted normal distribution based on the DFT analysis data. From the mathematical point of view, the normal distribution function is a form of probability function, hence, its cumulative expression takes the values from zero to unity. With this consideration, the $\alpha$ parameter is then the value of directly condensed adsorbate due to the surface curvature as explained by the Kelvin Equation (16). Therefore, if we compare the $\alpha$ values to the mesopore volume of the individual samples, a simple linear dependence is achieved in Figure 9, where the linear regression equals to $\alpha =0.8028d-0.0956$. This result can also be considered as a secondary proof of the model applicability. The only sample which does not fit to the equation is the SBA-15 sample, which was also excluded from the linear regression. The explanation could be that the different character of the SBA-15 adsorption makes the $\alpha$ parameter more dependent on the chemical properties of the material such as inverse influence on the number of layers described by the n constant. The SBA-15 has the highest layer number of $n=13.81$, while the other surface characteristics are comparable to the other samples. The SBA-15 also has significantly bigger pore size allowing for more surface coverage before the capillary condensation is in effect. However, to further study this dependence, more samples with larger pore size should be examined in a separate work.

The position of the capillary condensation on the x axis of the relative pressure $P/P_0$ is decided in Equation (33) by the $\omega$ constant. Furthermore, it is clear that the capillary condensation is distinguishably decided by the pore diameter. Therefore, considering the fact that from the definition of the Kelvin equation the critical radius, or curvature, is practically dependent on the natural logarithm of the relative pressure, reciprocally, the $\omega$, as a function of the relative pressure, will become dependent on the pore radius or diameter in this case. The linearisation through the logarithm function ($\omega =0.399ln\left(d\right)+0.0577$) shows a good fit with the data in Figure 10 for different pore diameters almost regardless the surface characteristics. This is, in fact, surprising as we have expected stronger dependence on the layer thickness, which directly affects the effective critical radius during the capillary condensation. Looking at the adsorption data presented in Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 it is apparent that the samples with aluminium have significantly higher adsorption uptake prior to the capillary condensation. A reasonable explanation can be provided by the mechanism of the adsorption on the silica surfaces enhanced for higher hydrophilicity by Brönsted acid sites such as the aluminium atoms incorporated in the silica structure. Even though that these atoms increase the adsorption potential of the silica materials, the effect is strongly localised and typical by cluster formations [63]. Based on this mechanism, the clusters despite significantly increasing the initial adsorption quantity have limited effect on the entire layer thickness, which consequently allows for a mostly undisturbed dependency of the $\omega$ on the real pore diameter. However, this clustering effect has also other implications which are moreover connected to the $\sigma$ constant.

The $\sigma$ parameter describes the span of pore capillary condensation, in other words, the range of the effective pore diameters associated with the capillary condensation. From the mathematical meaning, the smaller the $\sigma$ is, the narrower the effective pore diameter range will be. Accordingly, considering the clustering effect, which in terms of position is mostly random, the creation of localised water molecular bridges can be expected, broadening the capillary condensation part and increasing the $\sigma$ value. A similar effect of clustering can be observed in the case of microporous space (<$2\mathrm{n}\mathrm{m}$) which moves the equilibrium of adsorption dynamics further from desorption. This is clearly apparent from Figure 11, where the $\sigma$ value is shown as a function a natural logarithm of the ratio of C over n ($\sigma =0.0122ln\left(C/n\right)+0.0191$). Although, the selection is somehow limited, the materials with aluminium content or in the range of microporous region show notably higher $\sigma$ values. Furthermore considering the cluster formation as a basic process during the early stages of the adsorption, the final layer number n is not necessarily actual layer count, but rather an approximation of the clustering effect. Then, because the surface characteristics are simultaneously described by the C and n constants the function for the $\sigma$ parameter has to take in account both of these constants. From the definition it is clear that the higher is the basic surface affinity expressed by the constant C, as in the case if TMPS-1.5A, the higher the results are going to be, considering low n value. It is therefore apparent that the parameter $\sigma$ is dependent on two main factors. One is inherently connected to the porous structure, which is defined as the pore size distribution. The second factor is then the surface characteristics and the adsorption process defined independently by the BET constants and the apparent effect of cluster-like adsorption.

6. Conclusions

The problematic of mesopore filling is very complex and depends on many variables in general, and it is, therefore, impossible to describe it by a single or simple formula. It is apparent that the fitting of the mBET+G model is not entirely ideal and the thermodynamic properties of the system for the capillary condensation part cannot be described by this method yet. However, the actual fitting shows outstanding results and accuracy, and it can be surely considered for use in simulations of adsorption systems. The new model provides sufficient accuracy together with reasonable simplicity through the combination of all three main adsorption process types and points to an exciting direction for further development.

The advantage of the proposed mBET model in comparison to the recent models can be seen in adding a physically meaningful part describing capillary condensation region which has been so far difficult to address directly. The proposed model follows fully all the main modes of adsorption in a simple and clear way following the basic principles of both adsorbate–adsorbent and adsorbate–adsorbate interactions. Furthermore, the significance of the separate capillary condensation mode in the new model brings the possibility to apply standard models at high relative pressures, where they were usually failing, and obtain additional meaningful information about the modelled adsorption system in a simple and clear way compared to the other models.

We have explained that the new mBET+G model does not necessarily perform the best with all the samples; however, the overall performance is clearly superior and more universal to the other models. Hence, it should serve as a better tool to compare broader variety of materials. Even though the performance has been tested only on the materials with high level of structural organization, we believe that many other porous materials can be fitted with this model with good results and precision. Moreover, due to the parallel complementarity of the individual parts of the mBET+G (as well as GAB+G), the models can sustain the physical meaning of the constants from the original theories and create a unique system that can hopefully find its use in many areas of adsorption and porous material analysis.