How Reproducible are Surface Areas Calculated from the BET Equation?

Porosity and surface area analysis play a prominent role in modern materials science. At the heart of this sits the Brunauer–Emmett–Teller (BET) theory, which has been a remarkably successful contribution to the field of materials science. The BET method was developed in the 1930s for open surfaces but is now the most widely used metric for the estimation of surface areas of micro‐ and mesoporous materials. Despite its widespread use, the calculation of BET surface areas causes a spread in reported areas, resulting in reproducibility problems in both academia and industry. To prove this, for this analysis, 18 already‐measured raw adsorption isotherms were provided to sixty‐one labs, who were asked to calculate the corresponding BET areas. This round‐robin exercise resulted in a wide range of values. Here, the reproducibility of BET area determination from identical isotherms is demonstrated to be a largely ignored issue, raising critical concerns over the reliability of reported BET areas. To solve this major issue, a new computational approach to accurately and systematically determine the BET area of nanoporous materials is developed. The software, called “BET surface identification” (BETSI), expands on the well‐known Rouquerol criteria and makes an unambiguous BET area assignment possible.


Main
The Brunauer-Emmett-Teller (BET) equation is arguably one of the most used equations in physical chemistry and porosimetry. Since its conception in the 1930s [1] to estimate open surfaces whilst working with non-microporous adsorbents of the time such as Fe/Cu catalysts, silica gel, and charcoal, it found widespread use in the characterization of synthetic zeolites. [2] Furthermore, it has gained considerable momentum following the discovery of more complex porous materials such as mesoporous silicas, [3][4][5] porous coordination polymers (PCPs), [6] metal-organic frameworks (MOFs) [7] and covalent organic frameworks (COFs). [8] Novel porous materials are of significant academic and industrial interest due to their applications in gas storage and separation, [9][10][11][12] catalysis, [13] sensing, [14,15] and drug delivery. [16] To assess their adsorptive properties, Langmuir was the first to relate gas adsorption isotherms to surface areas, assuming only monolayer adsorption. [17] This was in contrast to Dubinin's proposition of micropore volumes for microporous materials. [18] Langmuir's adsorption theory was later extended to multilayer adsorption, resulting in the titular BET model. Even though the BET theory was not developed for describing adsorption in the microporosity, the BET area is now the de facto standard for the characterization of any porous material. Indeed, it has been recognized by the International Union of Pure and Applied Chemistry (IUPAC) as "the most widely used procedure for evaluating the surface area of porous and finely divided materials." [19,20] Furthermore, it has been an International Organization for Standardization (ISO) standard for surface area determination since 1995. [21] This makes it, arguably, the most important figure of merit for porous materials, including microporous ones. Looking at the literature, it is clear that the idea of monolayer coverage or even the concept of surface area are necessarily idealized and therefore could be inaccurate descriptions for microporous materials. [22] Indeed, Porosity and surface area analysis play a prominent role in modern materials science. At the heart of this sits the Brunauer-Emmett-Teller (BET) theory, which has been a remarkably successful contribution to the field of materials science. The BET method was developed in the 1930s for open surfaces but is now the most widely used metric for the estimation of surface areas of microand mesoporous materials. Despite its widespread use, the calculation of BET surface areas causes a spread in reported areas, resulting in reproducibility problems in both academia and industry. To prove this, for this analysis, 18 already-measured raw adsorption isotherms were provided to sixty-one labs, who were asked to calculate the corresponding BET areas. This roundrobin exercise resulted in a wide range of values. Here, the reproducibility of BET area determination from identical isotherms is demonstrated to be a largely ignored issue, raising critical concerns over the reliability of reported BET areas. To solve this major issue, a new computational approach to accurately and systematically determine the BET area of nanoporous materials is developed. The software, called "BET surface identification" (BETSI), expands on the well-known Rouquerol criteria and makes an unambiguous BET area assignment possible.
IUPAC warns users to apply "extreme caution [when using the BET equation] in the presence of micropores. (…) [The BET area] represents an apparent surface area, which may be regarded as a useful adsorbent 'fingerprint'." [19] This more nuanced understanding of the BET area is mirrored in the writing of Rouquerol et al., "the meaning of the BET surface is (…) that it embraces the major part of the amount of adsorptive in energetic interaction with the surface." [22] Despite these cautionary words, the BET area remains a deeply engrained metric in the fields of physical chemistry and materials science. Given the broad use of the BET equation, it is not surprising to see that much has been written on the applicability and the accuracy of the BET theory-that is, its model of the adsorption process-and on the reproducibility of the raw data, i.e., the adsorption isotherm. [23][24][25][26] Since the development of the first porous materials, there has been a sharp rise in the design of highly ordered and structured porous materials ( Figure S1, Supporting Information). [27,28] The advent of materials with more complex pore networks and dynamic frameworks through material design strategies such as reticular chemistry has given rise to reported BET areas in excess of 8000 m 2 g −1 . [10,[29][30][31][32][33][34] Often, these modern materials have complex adsorption isotherms which are more problematic or ambiguous to fit to the BET model, e.g. several steps can occur due to different pore types and/or flexibility being present in the material. [35] In response, a new generation of porosimetry equipment with pressure transducers capable of recording high-resolution gas adsorption isotherms at ultralow pressures (<10 −5 mmHg) has been developed. However, reliance on manual calculations of surface areas using the BET method remains commonplace. In this context, "manual" refers to the judicious selection of a pressure range by a scientist, be it through a self-developed spreadsheet or commercial software. This raises the question of the reproducibility of BET calculations from the same isotherm. An adsorption isotherm with 150 points has more than 10 000 consecutive combinations of points, all of which are potentially correct fitting ranges and will return different BET areas. The answer to the question of which is the optimal fitting region is far from obvious, and the consequence of any irreproducibility or different interpretations are serious. Consider two groups synthesizing the same compound and reporting two different BET areas; Sample A is reported to have a BET area of 1500 m 2 g −1 and Sample B's reported BET area is 2000 m 2 g −1 . Unless there is a common standard and protocol for calculating BET areas, we cannot say for certain that the quality and adsorption performance of Sample A is lower than that of Sample B. Indeed, the lack of reproducibility of MOF syntheses and adsorption performance, by comparing reported BET areas, has been highlighted already, but the natural spread of BET calculations was not included in the analysis. [36] The eponymously named Rouquerol criteria (Section S2, Supporting Information) aim to ensure good practice in identifying a valid fitting range, and, as such, they have found widespread acceptance in the literature and have been adopted in both IUPAC and ISO standards. [19][20][21][22]24,25,37] Despite this safeguard, we herein propose that current BET area calculations are irreproducible for two reasons: first, the Rouquerol criteria are indeterminate in identifying the correct fitting region, as they apply to multiple regions simultaneously. Second, even if they were determinate, they are too cumbersome and lengthy to implement and are therefore often neglected in practice. This dilemma is reminiscent of the Skeptic's Argument from Gorgias, here paraphrased: i) the BET area does not exist (e.g., for microporous materials); ii) even if it exists, it cannot be systematically and unequivocally calculated (i.e., determined by the Rouquerol criteria).
To prove our hypothesis and to assess the current spread of BET calculation results, we have shared a dataset of 18 isotherms (reported as relative pressure versus amount adsorbed), already measured, and representing four classes of micro-and mesoporous materials (zeolites, mesoporous silicas, MOFs, and COFs) with 61 laboratories with expertise in adsorption science and synthesis of porous materials. In this round-robin exercise, we asked the researchers to calculate the BET areas in the way they saw most fit. More details about the specific materials and the adsorption isotherms, sampled both from our laboratory and from the NIST/ARPA-E database, [38] are included in the Supporting Information, (Section S1). To avoid any recognition bias, all isotherms were anonymized and scaled-off arbitrarily. Figure 1a shows the large spread of results obtained from manual calculations of BET areas in the round-robin experiment, the full details can be found anonymized in the Supporting Information: tabulated in Section S3, and represented graphically in Section S4. Most groups (90%) reported using the Rouquerol criteria in their manual calculation, 23% used a commercial software package, and 6% used a self-developed code. Details on the methods for each group can be found in Section S11 of the Supporting Information. Bar a few exceptions, virtually no two groups of experts reported identical BET areas for any given isotherm. We observed a spread of at least 300 m 2 g −1 for each isotherm; however, that number was significantly higher for some individual isotherms. For NU-1104, a modern MOF with substantial porosity (isotherm Figure 1b), [32] the highest estimate of 9341 m 2 g −1 and the lowest estimate of  1757 m 2 g −1 differed by an astonishing 7584 m 2 g −1 , making the highest estimate more than five times higher than the lowest estimate. The spread of values for frequently reproduced MOFs such as HKUST-1, MOF-5, and ZIF-8 was slightly smaller than that of literature cited values. [36] While this observation affirms the natural assumption that material synthesis and isotherm measurement play a more important role in determining the BET area than the calculation, we can nevertheless attribute a significant portion of this spread to current BET fittings. The results of this social study demonstrate that there is significant variation in BET area calculations from the same isotherm, as it is extremely unlikely for two researchers to select identical fitting regions. At this point, we propose a novel method that not only systematically selects an optimal fitting region but does so by eliminating all other hypothetical fitting regions. With thousands of consecutive combinations of points, the large number of potential fittings is impossible to carry out manually.
To solve the problem of manual BET fitting, we developed a computational tool for BET analysis, BET surface identification (BETSI). This tool makes an unambiguous calculation of the BET area based on the original Rouquerol criteria but modified to prevent manual interaction, requiring only the adsorption isotherm as input data. As such, the results obtained from the round-robin evaluation were compared with the BETSI calculations to assess the inter-rater reliability of manual BET calculations.  Working principle for the BETSI algorithm. a) The isotherm is shown with a particular fitting region highlighted in red. The linear BET equation is applied, and an ordinary least-squares regression is applied to the fitting region. b) Subsequent checks against the Rouquerol criteria [24] are performed (insets) and c) valid fits are passed. The analysis shown in (a) is repeated for all consecutive combinations of points on the isotherm. A results matrix with n × n dimensionality stores all acceptable and rejected fits. d) All acceptable BET areas are output and plotted against the percentage error under the 4th Rouquerol criterion ((a), top inset). All BET areas ending on the highest permissible point under the 1st Rouquerol criterion ((a), bottom inset, maximum in the N(1 − (P/P 0 )) function) are labelled as the isotherm knee and shown in blue. The BETSI Optimal BET area (yellow) belongs to the isotherm knee group and has the lowest percentage error under the 4th Rouquerol criterion. algorithm on a simplified N 2 adsorption isotherm at 77 K for ZIF-8 (full details can be found in the Supporting Information, Section S5). First, the linearized BET equation is fitted to a particular region of the isotherm using an ordinary least-squares (OLS) regression (Figure 2a). The top panel shows the isotherm with a fitting region highlighted in red, and the OLS regression is shown below. The plot insets show the checks against the Rouquerol criteria (Figure 2b). If all criteria are met, the fitting is passed. This calculation is looped over all data intervals of at least 10 points on the isotherm. The resulting BET fits are stored in a large n × n matrix, where the (j,i)-matrix element corresponds to a fitting region starting at the jth-point and ending on the ith-point (Figure 2c). All valid fitting results are output and plotted against the percentage error under the 4th Rouquerol criterion (Figure 2d). Alongside, BETSI outputs all other BET parameters, such as monolayer capacity and the C constant, as well as full regression diagnostics (Section S5, Supporting Information).
Since multiple fittings comply with the Rouquerol criteria (Figure 2c,d), BETSI demonstrates that an unambiguous assignment of the BET area is impossible under the Rouquerol criteria alone. This proves our hypothesis that the criteria in their current form are indeterminate. For the prototypical ZIF-8 isotherm, a flexible MOF with narrow windows, [35] valid BET areas fall within a range of 1550 and 1750 m 2 g −1 (Figure 2c,d). BETSI assigns special relevance to fitting ranges that end on the highest permissible point, which are usually dictated by the 1st Rouquerol criterion, and labels these as the isotherm knee.
Beyond the isotherm knee, adsorptive activity decreases rapidly as the pores are mostly filled and the internal surfaces are saturated. Within this subset of BET areas, the BETSI optimum is chosen as the one with the smallest percentage error under the 4th Rouquerol criterion, thus making the BET assignment unambiguous.
Next, we ran BETSI on the isotherms distributed in the round-robin experiment. In all cases, the spread of potential BETSI results (i.e., those in agreement with the Rouquerol criteria) was considerably narrower than that obtained by manual calculation ( Table 1). Figure 3a shows the individual results from the social experiment and the comparison with the BETSI results; the corresponding variation coefficients are shown in Section S6 (Supporting Information) and an alternative representation normalized to the BETSI range is shown in Section S7 (Supporting Information). Since most groups reported using the Rouquerol criteria to calculate their BET areas, this substantiates our second hypothesis-that the manual implementation of the Rouquerol criteria is cumbersome and difficult to carry out in practice. For instance, in the case of NU-1104, the range of estimates decreases from 7500 m 2 g −1 in the social study to 235 m 2 g −1 under BETSI. Interestingly, some isotherms gave much larger spreads of results than others, suggesting that the BET model does not describe them as naturally and thus they are more susceptible to problems associated with the Rouquerol criteria. Unsurprisingly, we also observed this trend in the round-robin evaluation. To further investigate the goodness of the isotherm fittings, we define the BETSI variation Adv. Mater. 2022, 34, 2201502 Table 1. Results of BETSI analysis and round-robin evaluation for the isotherms used in the study. Material, isotherm of material under investigation; BETSI, optimal BET area predicted by BETSI; BETSI range, full spread of BET areas that pass under BETSI; BETSI variation coefficient, relative standard deviation of BET areas that pass under BETSI; pass rate, number of BET areas that pass under BETSI expressed as a fraction of all hypothetical fittings; round-robin average, mean of BET areas calculated in round-robin evaluation (sample size for each material, n = 61); round-robin range, full spread of BET areas determined in round-robin evaluation; round-robin variation coefficient, relative standard deviation of BET areas calculated in round-robin evaluation; hit rate, fraction of BET areas calculated in the round-robin evaluation that lie within the BETSI range.  Figure 3b demonstrates the correlation between the pass rate, the BETSI variation coefficient, and the hit rate. Simply put, the more BET fits are valid, the greater the spread of possible BET areas is, and the more likely researchers are to satisfy the Rouquerol criteria in manual calculations. To account for the non-equal spacing of points on all different isotherms, the pressure-adjusted pass rate expresses the total sum of pressure intervals that fit Rouquerol criteria as a fraction of the sum of all pressure intervals of the hypothetical fitting ranges (Section S8, Supporting Information). From Figure 3b, we classify adsorption isotherms into three broad categories, types A, B, and C (Figure 3c). While it is difficult to generalize about the shape of these isotherms, we still offer some discussion of common features.  Social study results versus BETSI results. a) Distribution of BET areas for identical isotherms from the social study (red) and BETSI (blue) obtained by kernel density estimation (sample size for each material, n = 61). Superimposed is the BETSI optimum (yellow). Note that the distributions of values obtained by BETSI are considerably narrower in all cases than those in the social study. b) Plot of the BETSI variation coefficient (relative standard deviation of BETSI results) against the pass rate (fraction of valid fits against all hypothetical ones). Bubble size scales with the hit rate, the fraction of results from the social study that lie within the BETSI range. The red symbols have a hit rate of zero. Note the positive correlation between all three parameters. c) Isotherm fit classifications. Type A fits have a relatively wide fitting window, within which multiple fits are possible, but return a relatively narrow spread of BET results. Type B fits have a narrow fitting window and concomitantly return a narrow set of spread of results. Type C fits have wide fitting windows, which translates to multiple passable fits and a wide spread of permissible BET areas.
computational support. In contrast to type A isotherms, type B isotherms often have sharp isotherm knees following strong adsorptive interactions at low relative pressures. Isotherms with more complex shapes such as NU-1104 also appear in this category. Type C isotherm fittings are arguably the most problematic. They have high pass rates and, concomitantly, they return large spreads of BET results. Typical materials that fit into this category are MIL-101, MIL-100, TPB-DMTP-COF and PCN-777. Like type A isotherms, these have rounded isotherm knees, which appear at higher relative pressures. It is for these materials that the necessity to extend the Rouquerol criteria is demonstrated and the BETSI algorithm makes an unambiguous BET assignment possible.

Outlook
BET theory is a great success story. Developed in the 1930s for non-microporous, open surfaces, it continues to this day to be applied to modern adsorbents with complex porosity. Despite the advances from classical density functional theory (DFT) methods, the BET area will likely continue playing a crucial role in porosimetry for decades to come, with impacts in energy research, transport, medical applications and climate-change mitigation. In light of these future developments, it will become increasingly important to share critical scientific metrics reliably to find a common language to report both academic and industrial progress.
Here, we have demonstrated the difficulties in unambiguously determining BET areas from adsorption isotherms, which, in turn, affect the assessment of material quality and reproducibility. These problems arise from imperfect and insufficient manual calculations and can only be met using modern computational methods. Furthermore, we propose BETSI as a step toward greater transparency and criticality in determining BET areas. We stress here that it is neither the function nor the purpose of BETSI to eliminate doubt and treat a particular BET area as "true." Researchers should remain aware of the limitations of BET theory when applied to microporous adsorbents in general and when BET areas are reported, the pressure range and number of points used should always be stated. We further recommend here that isotherms must be reported transparently and in detail, i.e., semi-log representation to show the low-pressure regions. The "experiment" is the adsorption isotherm-not the BET area.

Round-Robin Evaluation
N 2 adsorption isotherms of 18 different materials (Section S10, Supporting Information) were sent to international collaborators: HKUST-1, ZIF-8, NU-1000, MIL-101, UiO-66, Al fumarate, Zeolite13X, Mg-MOF-74, UiO-66-NH 2 , MOF-5, DMOF-1, MCM-41, TPB-DMTP-COF, MIL-100, NU-1102, NU-1104, NU-1105, and PCN-777; they were anonymized and labelled A-R respectively. Note that this is not the order in which the isotherms appear in the paper. The isotherms were sampled from the A 2 ML measurements and from the NIST Adsorption Database. Arbitrary scaling factors were introduced to minimize recollection bias of the isotherms. The isotherms were sent out in .csv format. All colleagues received the same email with the same set of instructions (Section S1, Supporting Information): to calculate the BET area from the data in the way they saw most fit and to report a rough estimate of how long it took them to calculate them. An anonymized one-page summary of each lab's own account of their calculation can be found in Section S11 (Supporting Information).
For easier data handling, once rescaled, all results were rounded to the next integer. None of the data points has been eliminated. The data is presented as a jitter plot for each material, with a superimposed kernel-density estimation obtained in Python.

BETSI
The BETSI algorithm, including executables (for Linux, Windows, and macOS) and a user manual are available at ref. [39]. The "How to Use BETSI" guide is included in the Supporting Infomation (Section S15). The program is written in Python and uses principally the NumPy library. It loops linear regressions over all consecutive combinations of at least three points, performs full BET analyses, and stores the fitting parameters in n × n results matrices, where the (j,i)-matrix element denotes a linear regression from the jth to the ith point on the isotherm. Binary pass/fail matrices with the same dimensionality are used independently to assess compliance with linearity and fitting criteria. The "filtering" of BET areas is achieved by elementwise matrix multiplication of the results matrices and the pass/ fail matrices. This allows independent "activation" and "deactivation" of the criteria and observing the effects on the results. The minimum fitting requirement of ten points is coded in a pass/fail matrix to allow for some minimum point flexibility, as is the cut-off value for R 2 of 0.995. To avoid low-leverage nonlinearity in the linear region, the first Rouquerol criterion has been extended to also require the linearized BET function to increase monotonically with P/P 0 , as well as N(1 − (P/P 0 )). The third and fourth Rouquerol criteria are implemented through a 10 000-point Pchip interpolation of the isotherm to reconstruct the N m (Read). As the third and fourth criteria require the N m (BET) to be a real value, i.e., they require C to be positive, the second criterion cannot be independently deactivated from the third and the fourth. The associated logic has been written into the program. Following the BETSI filtering by multiplication of results and pass/fail matrices, the isotherm knee is identified as the subset of BET areas whose fitting region end on the highest permissible pressure point. In most cases, this will be the highest permissible point under the first Rouquerol criterion. The optimal BETSI prediction is chosen as the fitting region with the lowest percentage error under the fourth criterion and belonging to the isotherm knee subset.
BETSI only requires the adsorption isotherm as input data and returns six plots used to validate the results: the isotherm itself, with the optimal linear region highlighted as well as the BET fit; the "Rouquerol representation" of the isotherm, N(1 − (P/P 0 )) plotted against P/P 0 ; the linearized plot with the OLS regression and the regression parameters; the filtered percentage error versus BET areas plot with the isotherm knee and optimal BET area highlighted; the filtered monolayerloadings plot showing all permissible monolayer loadings on the isotherm; and the statistical distribution of permissible BET areas with a boxplot. Additionally, BETSI returns four regression diagnostics plots, which can be used to assess whether the assumptions of OLS regression have been met: the residuals versus fitted values plot can be used to visually inspect whether the residuals are normally distributed around the regression line, and similar information could be obtained from the QQ-plot. Finally, the scale-location plot can be used to assess whether the distribution of studentized residuals is homoscedatic or heteroscedatic and the residuals versus leverage plot can be used to identify high-leverage points that have an abnormally large influence on the regression line.

Comparison between Round-Robin Evaluation and BETSI Results
Statistical analysis of the results was performed in Python. The BETSI variation coefficient and the round-robin variation co efficient are standard deviations relative to the average of each set. The pass rate for each isotherm is the number of permissible BET fits as a fraction of all consecutive combination of points. To account for non-equal spacing of the points on each isotherm, the pressure-adjusted pass-rate is obtained by integrating along the pressure axis and dividing the total sum of permissive pressure intervals by the sum of all consecutive pressure intervals. The hit rate is the fractional number of BET areas calculated in the round-robin evaluation that lie within the BETSI range.

Statistical Analysis
Data was analyzed using the NumPy, SciPy, and Matplotlib libraries in Python. For the violin plots, a sample size, n = 61 was used for each material, and the estimator bandwidth for the kernel density estimator (KDE) was calculated using Scott's rule. For the bubble plots, the hit rate has been averaged over a sample size, n = 61, for each material.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.