Water vapor adsorption-desorption hysteresis due to clustering of water on nonporous surfaces

Water vapor adsorption-desorption hysteresis due to clustering of water on nonporous surfaces Yrjö Viisanen1, Maria Lbadoui-Darvas2,3,4, Ana A. Piedehierro1, André Welti1, Athanasios Nenes2,3, Ari Laaksonen1,5* 1Finnish Meteorological Institute, 00560 Helsinki, Finland. 2Laboratory of Atmospheric Processes and their Impacts, School of Architecture, Civil & Environmental Engineering, École Polytechnique Fédérale de Lausanne, Lausanne, 1015, Switzerland 3 Institute of Chemical Engineering Sciences, Foundation for Research and Technology (FORTH/ICE-HT), Hellas, 26504, Patras, Greece 4Federal Office of Meteorology and Climatology MeteoSwiss, Chemin de l'Aerologie 1, CH-1530 Payerne, Switzerland 5Department of Technical Physics, University of Eastern Finland, 70211 Kuopio, Finland. *Correspondence to: [email protected] Abstract Water vapor is continuously adsorbed onto and desorbed from all kinds of surfaces depending on changes of relative humidity. Adsorption-desorption hysteresis of water that occurs on various nonporous surfaces and extends down to low relative humidities has been reported for decades but remains unexplained. Here we show experimentally that such hysteresis is a common phenomenon on metal oxide and mineral surfaces and can be divided into two distinct categories based on the wettability of the adsorbent surface. Type I hysteresis occurs on more hydrophobic surfaces and is associated with adsorption isotherms that behave rather linearly as water saturation is approached, whereas type II hysteresis occurs on more hydrophilic surfaces and is associated with adsorption isotherms that curve steeply upward close to saturation. Our model calculations provide strong indication that adsorption in both types occurs cluster-wise, and the type I hysteresis is caused by contact angle hysteresis while type II hysteresis is associated with film formation close to saturation. The understanding of water vapor adsorption and desorption mechanisms may be key for explaining and quantifying physical and chemical interfacial phenomena in atmospheric and industrial environments. INTRODUCTION Water vapor is readily adsorbed on practically all surfaces at ambient conditions. Depending on relative humidity (RH), temperature, and chemical and physical properties of the surface, adsorbed water can exist as individual molecules, molecular clusters, or monoor multilayer liquid or ice-like films (1). Interfacial water impacts many industrial and environmental processes, including atmospheric chemistry. In this context, water vapor adsorption on both hydrophilic and hydrophobic surfaces has been extensively studied for more than a century. It is therefore rather surprising that the cause of a phenomenon first observed already decades ago (2-9), the adsorptiondesorption hysteresis of water extending down to low relative humidities, has so far remained poorly understood. Adsorption-desorption hysteresis is a commonplace phenomenon in porous materials (10). Capillary condensation, i.e. filling of pores by the adsorbate, occurs at vapor pressure lower than the saturation vapor pressure of the adsorbate due to concavely curved surfaces of the pores. Evaporation from the filled pores does not take place reversibly but requires a lower saturation ratio S than for condensation, hence a hysteresis loop is seen in a plot of adsorbed amount vs. saturation ratio (with water vapor, S = RH/100%). The capillary condensation and evaporation saturation ratios depend on pore sizes and geometries; however, the hysteresis loop is always closed, i.e. the hysteresis does not extend to very low saturation ratios. For example, with N2, the hysteresis does not extend to saturation ratios below 0.4 – 0.5 The exact saturation ratio at which the loop closes is dictated by the sizes of the smallest pores in the adsorbent material into which the adsorbate molecules can penetrate. Fig. 1 shows an example of water vapor adsorption hysteresis in mesoporous silica; in this case the hysteresis loop closes at 30% RH, i.e. water vapor saturation ratio of 0.3. Figure 1. Experimental water adsorption and desorption isotherms showing hysteresis between 30-70 % RH. Apart from pores, adsorption-desorption hysteresis can also result from the effects of water chemisorption (11). Water molecules dissociate to OHand H+ on many surfaces (1) and bind more strongly than water molecules themselves. For example, on silica surfaces this so-called hydroxylation process leads to formation of silanol (SiOH) groups (12, 13) that bind water molecules effectively so that a hydroxylated silica surface can adsorb more water vapor via physisorption than a bare silica surface (14). The silanol groups can only be removed by heating the surfaces to sufficiently high temperature, whereas physisorbed water molecules can be removed from the surfaces at room temperature by vacuum pumping. Nevertheless, chemisorption-caused hysteresis will not be observed once the surface has been saturated with respect to chemisorption, or fully hydroxylated (15), i.e., after an adsorption-desorption cycle (11) during which the hydroxyl formation takes place on all possible surface sites where it can occur. With environmental processes, the adsorbent surfaces are often exposed to many cycles of humidification and dehumidification, so chemisorption is not expected to have impact on the related adsorption and desorption curves. From the above, the water adsorption-desorption hysteresis loops that extend to very low saturation ratios (after the surface has already gone through one or more adsorption-desorption cycles) are neither caused by capillary effects nor by dissociation of water molecules. Clustering of adsorbed water (1, 16) however has been suggested to cause the hysteresis (8). Here, we provide experimental evidence that the hysteresis takes place on many different nonporous surfaces and show, using molecular simulations and theoretical considerations, that clustering of water molecules can indeed lead to the observed differences between the adsorption and desorption isotherms. Thus, an inherent, ubiquitous but ignored adsorption mechanism – that does not require pores or chemical interactions with a surface – may on its own provide the interactions required for explaining and predicting aerosol-water interactions in atmospheric and other applications. We furthermore suggest – supported from experimental evidence – that there are in fact two types of hysteresis loops, type I occurring on less wettable and type II on more wettable surfaces. We propose that the type I hysteresis is caused by the difference between advancing and receding contact angles of water clusters during adsorption and desorption, respectively, whereas type II hysteresis is caused by water clusters fusing into a film along the adsorption branch. Our type I and type II hysteresis are not to be confused with IUPAC classification of hysteresis loops H1-H5 caused by different types of pores and pore networks. EXPERIMENTS AND MODELLING Experimental Water adsorption isotherms were measured volumetrically at room temperature using a Belsorp MAX II automatic gas-adsorption instrument (Microtrac MRB). Table 1 shows the materials used in the adsorption measurements. The specific surface areas of the samples were determined via BET analysis (17) of nitrogen adsorption isotherms measured using the same setup. The statistical number of adsorbed water monolayers N was calculated based on the N2 surface areas, assuming adsorbed nitrogen and water molecular cross-sectional areas of 16.2 and 10.5 Å2, respectively (18). The absence of hysteresis in the nitrogen adsorption isotherm measurements (see SI) indicates that the analyzed samples are non-porous. Several water adsorption isotherms for the same sample are available (not included here) showing the repeatability of the measurements. Additionally, prior to the measurements, a series of water adsorption cycles were performed with every sample, ensuring that chemisorption of water molecules reached a saturated state. This guarantees physisorption as the only adsorption mechanism taking place on the samples during the subsequent measurements analyzed in the present paper. Table 1. Summarized information about the powder materials used in this study. Grain size is provided when available. The surface areas are derived from N2 adsorption measurements. Isotherm type Material Producer Purity Surface area (m2g-1) I Aluminum oxide < 1.0 micron powder Al2O3 Alfa Aesar 99.9% (metal basis) 4,48 I Copper(II) oxide CuO Alfa Aesar 99.995% (metal basis) 0.27 I Iron(III) oxide Fe2O3 Alfa Aesar 99.9% (metal basis) 2.30 II Kaolin Al2Si2O5(OH)4 Riedel-de Haën Natural sample 9.12 II Nickel(II) oxide 325 mesh powder NiO Alfa Aesar 99% (metal basis) 0.54 I Lead(II) oxide PbO Alfa Aesar 99.99% (metal basis) 0.55 II Lead(IV) oxide PbO2 Alfa Aesar 87.0% (min) 1.43 II Lead(II,IV) oxide Pb3O4 Alfa Aesar 97% (metal basis) 0.67 II Tin(IV) oxide < 10 micron powder SnO2 Alfa Aesar 99.9% (metal basis) 6.36 I Titanium(IV) oxide 0.9-1.6 micron ASP powder TiO2 Alfa Aesar 99.8% (metal basis) (min 97% rutile) 4.41 I Titanium(II) oxide 325 mesh powder Ti2O3 Alfa Aesar 99.8% (metal basis) 0.98 Molecular simulations Sessile droplet molecular dynamics (MD) simulations on graphitic slabs were performed using GROMACS 2022 (19) to study the contact angle hysteresis of water clusters consisting maximally of some thousands of molecules. The carbon atoms in the graphitic slab were kept immobile during the simulations by means of position restraints with force constants of 1000 kJ/mol/nm in each direction. Water molecules were described using the SPC/E (Simplified point charge/Extended) (20) water model and graphite-water interactions were modelled using the Lennard-Jones (LJ) parameter set proposed by Werder et. al. (21) that maximizes the agreement between simulated and experimental contact angles (ε=0.3920 kJmol-1 ς=0.315 nm). The Lennard-Jones interaction potential between the centers of two atoms at distance (r) can be written as (22): 𝑉𝐿𝐽(𝑟)=4𝜀[(𝜍 𝑟)12 − (𝜍 𝑟)6]. The geometry of the water molecules was kept constant using the LINCS algorithm (23). LJ interactions were truncated to 0 beyond an atom-based spherical cutoff of 1.2 nm and long range Coulombic inter-actions were approximated by the Particle Mesh Ewald (PME) algorithm (24) beyond the same cut-off distance. All simulations were performed on the canonical (NVT) ensemble. The temperature was controlled using the V-rescale thermostat (25). Each simulation was initiated from a levitated pre-equilibrated spherical water droplet over a graphite slab consisting of 10 layers with dimensions of 25 by 25 nm. In the equilibration phase, common to all simulations, a water droplet consisting of 1000, 7000 or 20000 molecules was allowed to settle on the surface in a 10 ns long equilibrium simulation on the canonical (NVT) ensemble. This method ensures that the droplet contact angle is not biased by the initial configuration being in a metastable state, which can occur if the droplet is posed directly in contact with the surface at a pre-set contact angle (26). The length of the equilibration simulation was determined by monitoring the changes in the contact area between the droplet and the surface. The system was considered equilibrated when the contact area reached a plateau. In order to monitor the dynamic changes in the contact angle during evaporation, simulated annealing was performed during which the system was gradually heated at an annealing rate of 0.4 K/ns. The initial temperature was 298 K and the final was 500 K. These simulations were performed only with the largest droplet containing 20000 water molecules. In order to simulate the quasi-equilibrium evaporation, model equilibrium MD simulations were performed using the equilibrated systems with the three different droplet sizes at T=298 K. These latter simulations were 50 ns long, the total of which was used to calculate the contact angle. We note that this is a standard process to study contact angles of fluids on surfaces (26), however it does not allow setting the relative humidity in an explicit manner, i.e., by adding a fixed number of water molecules to the gas phase. The fast diffusion of and the large distances between molecules in the gas phase create a situation in which zero water molecules in the gas phase represented by a limited empty volume in our simulation box can correspond to the same RH as multiple water molecules. Rephrased in the framework of statistical thermodynamics, these configurations with zero or multiple gas phase molecules are different microstates in the ensemble described by a given temperature, total number of water molecules and pressure (or RH). Therefore, this simulation corresponds to any RH that permits the formation of a water droplet on the surface of graphite. The time evolution of the average contact angle was calculated as described below from the height and the radius of the spherical cap for the annealing simulations, however this method is only valid for larger droplets (R > 5 nm) where the droplet assumes a spherical cap geometry and perturbations from the fluctuations of the droplet surface are negligible compared to the mean droplet radius (26,27). For the equilibrium MD simulations, a more exact method was developed which relies on distributions of local contact angles in small elements of the droplet surface (see below). The equilibrium contact angles were estimated as the median of the obtained contact angle distribution. Calculation of the equilibrium contact angle for the annealing simulations The time evolution of the contact angle was calculated from the annealing simulations using the base radius (b) and height (h) of the droplet in the spherical cap approximation (28). In this approximation the shape of the droplet was estimated as a spherical cap whose height (h) is the vertical distance between the position of the plane of the top graphite layer and the water molecule furthest from it in the Z direction (which is vertical to the surface). The base radius (b) of the spherical cap was taken as half of the mean maximal distance between the water molecules constituting the contact layer with graphite slab in the X and Y directions (which are parallel to the surface). In this framework the contact angle (Θ) can be calculated from cos Θ = (R-h)/R where R=h/2+b2/2h. Figure 2. shows the definition of the variables needed for the calculation. Figure 2. Schematic representation of the data used to calculate the contact angle of a spherical cap. Note that for droplets whose radius is smaller than 5 nm, this approximation is crude because the fluctuations of the liquid surface are comparable with the droplet radius, whose shape therefore deviates significantly from a spherical cap. Calculation of the contact angle from the local contact angle distribution We developed a method to estimate the local droplet shape based on the triangulation of the droplet using the surface molecules as apices. It relies on selecting the molecules that constitute the surface at each frame of the simulation using the GITIM algorithm (29) implemented in the PYTIM library (30). Molecules that constitute solid/liquid surface – i.e. those that are within 0.35 nm from the graphite surface – are also identified as interfacial by the GITIM algorithm, and are excluded from further analysis. Once the surface molecules are selected and the solid/liquid surface is excluded, we perform the triangulation by identifying neighbor triplets using the heapsort (31) algorithm to simulation. The experimental contact angle decreases almost linearly from slightly above 80° degrees to slightly above 60° by the time the droplet volume has decreased by 40%. Thereafter, the decrease is less steep, and the final value for a droplet with 10% water left is slightly under 60°. The first simulation method (simulated annealing) models the dynamic change in the droplet volume because of evaporation. The simulation was initiated with a droplet whose original radius was 7 nm at T=298 K. As a result of the gradual heating, the droplet began to evaporate and shrink. The contact angle of the droplet in the beginning of annealing phase was approximately 95°. At the end of the simulation, when the droplet shrank to about 75% of its original volume, the contact angle was approximately 75°. The rate of contact angle change is higher in the simulations than in the experiments between 90-100% droplet volumes (i.e. in the beginning of the evaporation), then the simulated curve changes slope at lower droplet volume percentages and aligns with the experimental curve. The deviation in the initial phase potentially results from the fact that annealing simulations - due to limitations in computing power - were performed at a very high heating rate (0.4 K/ns), thus the resulting evaporation is strongly out of equilibrium. At a few nanoseconds temporal resolution, evaporation due to vacuum pumping can be considered as a series of quasi equilibrium steps. Equilibrium molecular dynamics simulations with three different droplet volumes representing 100%, 37% and 9% of the original 7 nm droplet were performed at constant temperature (298 K) to model the contact angle changes in the quasiequilibrium approximation. The equilibrium contact angle was calculated in these simulations using a local curvature scan algorithm (see above). The contact angles, 82° at 100% and 57° at 9% volume, agree well with the experimental curve while the one at 37% droplet volume is 10° higher than the experimental value. Overall, molecular simulations - both dynamic and quasi equilibrium approaches - confirm that evaporation results in a decreased contact angle even in nanodroplets that contain only a few thousands of molecules, which supports that the assumption of contact angle hysteresis being responsible for type I adsorption - desorption hysteresis is plausible. ab c d e f Figure 4. a) Equilibrium configuration of the molecular simulation of an evaporating water nanodroplet b) Receding contact angle of water on graphene as a function of droplet volume. c) Type I adsorptiondesorption hysteresis of water on model graphene. Red line: adsorption, grey line: desorption. The y-axis shows statistical number of water monolayers on graphene. d) Type II adsorption-desorption hysteresis of water on model silica. e) Type I adsorption-desorption hysteresis mechanism of water on graphene. The contact angle along the adsorption branch is higher than for the desorption branch. f) Type II adsorptiondesorption hysteresis mechanism of water on model silica, where a film forms at 92.2% RH. Desorption takes place from the film. The type 1 adsorption-desorption hysteresis loop for water droplets on graphite was modelled (Fig. 4c) using the same extended FHH adsorption parameters as in our earlier work (47), except that the experimental advancing and receding contact angle values were applied. Thus, the advancing angle along the adsorption branch was 82°; and the receding angle decreased with droplet size similarly as in the experiments (Fig. 4b), ending at a value of 57°. As expected, the adsorption isotherm touches the vertical axis at 100% RH, signalling that the growing droplets remain stable at saturation. In fact, the adsorption isotherm could be continued into the supersaturated region, until a maximum (“critical”) supersaturation is reached. This is a point where the droplets undergo heterogeneous nucleation (37) and drive the “unimpeded” condensation of water vapor onto the surface. The shape of the desorption isotherm differs from the adsorption isotherm in that, the radius of curvature of the desorption isotherm changes sign twice, at low RH below 5% and at around 60%; whereas the radius of curvature of the adsorption isotherm changes sign only at low RH. Fig. 4d shows type II hysteresis on a model surface. The adsorption curve (red) was calculated using the extended FHH model for hydrophobic silica. In practice, we used the same model parameters as in our earlier work (48) (A = 1.88; B = 1.54; Θ = 20°) but reduced the density of water adsorption sites by increasing the average distance between sites to 10 nm (instead of 1.61 nm used in (48)). Note that while this increases the hydrophobicity of the surface, capturing the real adsorption characteristics of hydrophobic silica with reduced density of silanol groups might require adjusting the FHH parameter and contact angle values as well as the adsorption site density. We however skip this exercise as our purpose here is not to reproduce the water uptake by real hydrophobic silica as closely as possible but to demonstrate qualitatively how the extended FHH model generates type II hysteresis. Droplet-wise adsorption takes place on the model silica surface up to a relative humidity of 92.2%, where the curve jumps up as the nanodroplets fuse and form a uniform film. In the real world, such stepwise features from film formation are likely smoothed out as clusters are expected to gradually fuse to larger patches before the complete film forms. The black desorption curve in Fig. 4d was computed using the traditional multilayer FHH model. Proceeding to our experimental results, the materials used in the experiments are presented in Table 1. A common feature of all water adsorption-desorption hysteresis loops is that they span practically the whole dynamic RH range (see Figures 5 a and b). Furthermore, the experimental isotherms can be divided into two types as suggested by our theoretical calculations. The adsorption isotherms of type I hysteresis loops (Fig. 5a) do not turn upward as 100% RH is approached, whereas those of type II loops (Fig. 5b) do. Note that the y-axes of the plots are logarithmic as the adsorbed volume of water on the different material samples varies by almost two orders of magnitude; individual isotherms are plotted in linear form in Figs. 5 c and d, and in figures S1 and S2 (see SI). It can be seen from Fig. S1 that the desorption isotherms of at least Fe2O3 and TiO2 resemble those of the modelled graphene (Fig. 4c) in that their radii of curvature change sign at intermediate RH (about 70%). Figures 5 c and d present comparisons between modelled and observed isotherms for both type I and type II hysteresis; the extended FHH model captures the observations quite well at RH above 20% (as noted above, at lower humidities, some of the model assumptions become overly simplified). For example, the water adsorption for PbO shows a wide hysteresis loop that could be reproduced using FHH parameters of A = 7.5 and B = 2.1 (Fig 5c). We found that applying an advancing contact angle of 92° with 2.25 nm average distance between nanodroplets resulted in a satisfactory adsorption isotherm, and a desorption isotherm matching the experimental data was obtained by letting the receding contact angle to decrease from 92° linearly to 82° as RH approached 0%. A phenomenon that may influence our experimental isotherms is capillary condensation at interparticle contact points within the powder samples. However, with powders consisting of super-micron particles (which is mostly the case with our samples), interparticle capillary condensation occurs at high relative humidities (RH > 98%) (49). Furthermore, if this phenomenon induced a hysteresis that extends to clearly lower relative humidities, it should be detectable in the N2 sorption measurements as well, but there is no hysteresis (see SI). The same arguments apply to possible hysteresis caused by percolation within the powders. However, it is possible that interparticle capillary condensation at high RH contributes to the film formation in type II hysteresis. Figure 5d shows a comparison of theoretical and experimental water adsorptiondesorption loops on kaolinite. The adsorption-desorption hysteresis loop appears narrower for kaolinite than for PbO, modelled graphene or silica, as is the case with some of the other adsorbents (see Fig. S2). The root cause for the narrow hysteresis loops is that in our adsorption experiments, we use powder samples, which in most cases contain more than one crystal faces of the adsorbent material. The different crystal faces are likely to have different water contact angles and FHH parameters. Molecular simulations of water on kaolinite (50) support this; kaolinite is a layered material composed of tetrahedral silicate sheets connected to octahedral alumina sheets via hydrogen bonds, and the simulations suggest that the octahedral surface is essentially completely wettable by water, whereas the tetrahedral surface is hydrophobic. The water adsorption isotherms have also been found to differ on the different crystal faces of kaolinite (51). For the extended FHH-theory calculations of adsorption and desorption of water on kaolinite, we assumed that 50% of the surface (octahedral) is completely wettable with film-wise adsorption, and the other 50% (tetrahedral) is hydrophobic with cluster-wise water adsorption. The FHH parameters A and B of the octahedral and tetrahedral surfaces (Ao, Bo and At, Bt, respectively) were selected so that the completely wettable octahedral surface takes up more water vapor than the hydrophobic tetrahedral surface (see SI for details of the FHH parameter selection). The selection of contact angle (86°) and average distance between water clusters (2.3 nm) at the tetrahedral surface that reproduced the adsorption isotherm well (see SI) resulted in a model where the clusters a b c d Figure 5. a) Experimental adsorption-desorption hysteresis curves of type 1. b) Experimental adsorptiondesorption hysteresis curves of type 2. c) Experimental and modelled hysteresis of type 1 for water adsorption-desorption on PbO d) Experimental and modelled hysteresis of type 2 for water adsorptiondesorption on kaolinite. on the adsorption branch, red curve in Fig. 5d, fuse to form a film at a relative humidity of over 99%. The film then persists throughout the desorption branch, grey curve in Fig. 5d, until a monolayer is left. It is evident that the modelled adsorption and desorption curves match the data well above about 20% RH. It should be noted that in our modelling, we have not included water adsorption on the edge surfaces of kaolinite sheets, which, according to the recent molecular simulations of Zhang and coworkers (52), can be significant especially at low relative humidities. Regarding adsorption on the octahedral and tetrahedral surfaces, our results appear consistent with the simulations (52). The good match of our model calculations with the experimental results provides compelling evidence of cluster-wise adsorption of water vapor. The physical properties of water clusters differ from those of water films, which may cause erroneous predictions by adsorption models in which film-wise adsorption is assumed. For example, we find a specific feature (explained below) in our extended FHH model calculations which we also see in our experimental data that casts doubt on the applicability of so-called BET surface area analysis using water adsorption isotherms. BET adsorption theory (17) predicts that at saturation ratios between 0.05 and 0.35, a plot of [V(S-1-1)]-1 vs S (where V is the volume of the adsorbed gas) is linear, and that the monolayer coverage Vm of adsorbate molecules on the surface is obtained from the slope and the yintercept of the plot. If the molecular cross-section of the adsorbate molecule is known, the specific surface area of the adsorbent can be calculated. BET analysis is mostly used for determining specific surface areas of powdered materials with nitrogen or some other inert gas as adsorbate; it is sometimes applied also for water isotherms (e.g. 53-59), and the commonly used value for the cross-section of water molecule, 12.5 Å2, results partly from BET analysis done for water isotherms on several adsorbents (18) (in contrast, the value we use in our model calculations, 10.5 Å2, comes from theoretical considerations). A BET plot calculated from the extended FHH results for the model silica of Fig. 4d is shown in Fig. 6a. It is seen that the resulting curve (black dots) is nonlinear (downward concave). This feature is produced by the curvature of the nanodroplets, with the extent of nonlinearity depending on the model parameter values. In contrast, the BET plot calculated from the traditional FHH isotherm (corresponding to the desorption isotherm of Fig. 4d), grey crosses in Fig. 6a, is almost perfectly linear. Fig. 6b shows the nonlinearity of the BET variable for the model silica and six of our experimental adsorbents. The BET variables of each adsorbent were divided by the corresponding values obtained from linear least square fits at relative humidities between 5 % and 35 %. The shapes of the resulting curves vary somewhat, but in all cases the downward concave feature predicted by the extended FHH theory is seen. a b Figure 6. a) BET plot for model silica. Black dots: data points calculated using the extended FHH model; red line: linear lest squares fit to the data points. Grey crosses: data points calculated using the traditional FHH model. The data points were shifted upward by 5 units to improve readability of the figure. Grey line: linear least squares fit to the data points. The extended FHH model values correspond to the adsorption curve and the traditional FHH model values to the desorption curve of Fig. 4d. b) Modified BET plots calculated from adsorption isotherms of water on six of the studied adsorbents and the model silica of panel a. The curves show the ratio (QBET) of the measured [V(S-1-1)]-1 to the respective linear least-square fit at relative humidities between 5 % and 35 %. The nonlinearity of the BET variable calculated from the experimental adsorption isotherm is, beside the hysteresis loop, an additional indication of cluster-wise adsorption. However, we do not observe the nonlinear feature with all our adsorbents. With adsorbents showing type 2 hysteresis, this may be due to adsorption at low humidities occurring mainly on wettable crystal faces. It is also possible that some combinations of adsorption parameters produce quite minimal nonlinearity. In any case, our extended FHH model calculations indicate that beside inducing nonlinear behavior of the BET variable, the results from the BET analysis in case of cluster-wise adsorption tend to lead to (sometimes strongly) underestimated specific surface areas. We therefore do not recommend doing BET surface area analysis using water adsorption isotherms. CONCLUSIONS Adsorption-desorption hysteresis of water vapor on various non-porous surfaces spanning the whole RH range has been reported by many researchers during the past several decades. However, no convincing explanation of the mechanism(s) producing the hysteresis has been offered. At least in two instances it has been suggested that clustering of water might cause the hysteresis (6,24), but without further evidence, the suggestions have remained speculative. In this work we used a combination of experimental and modelling approaches to demonstrate the very likely reasons for the adsorption-desorption hysteresis. We measured adsorption and desorption isotherms on a number of non-porous metal oxide and mineral surfaces, making sure that the surfaces were saturated with respect to chemisorption before the measurements. We found the hysteresis to be a common phenomenon and were able to group the measured isotherms into two categories based on the shapes of the isotherms. In type I hysteresis, the isotherms remain linear as 100% RH is approached, whereas in type II hysteresis, the isotherms turn steeply upward close to saturation. Our results suggest that the type I hysteresis is caused by contact angle hysteresis: in the adsorption branch, the contact angles of water clusters are larger than on the desorption branch, causing the adsorption branch clusters to be smaller at a fixed RH than the desorption branch clusters. Molecular dynamics simulations confirmed that the contact angle hysteresis is of similar magnitude for water clusters consisting of a few thousands of molecules (i.e. typical sizes of adsorbed clusters at high relative humidities) as for milliliter sized droplets. Making use of this information, we demonstrated with the extended FHH model that adsorption-desorption hysteresis similar as that seen in the experiments can be reproduced assuming reasonable contact angle hysteresis. We furthermore showed that the type II hysteresis is very likely caused by water clusters filling the complete surface somewhere along the adsorption branch. Desorption therefore takes place from a uniform film with zero curvature, whereby the statistical thickness of the water layer at a given RH (below that where the filling of the surface takes place on the adsorption branch) is larger than on the adsorption branch. We were able to model the type II hysteresis on kaolinite successfully assuming that the octahedral surfaces are completely wettable and the tetrahedral surfaces hydrophobic. Our extended FHH model predicts that in many cases, the BET plots become distinctly non-linear (downward concave) when adsorption is cluster-wise. A similar feature was seen in the experimental data with several adsorbents. The extended FHH model results suggest that the non-linearity can lead to underestimated specific surface areas if BET analysis is carried out. The quantitative modelling of the adsorption-desorption hysteresis of water vapor measured on powdered materials is somewhat difficult because powders tend to contain different crystal faces that may have different wetting properties and adsorption parameters. In this respect, adsorption and desorption isotherms measured on single crystal faces, or isotherms measured on amorphous (or glassy) powders, would likely be easier to characterize. In any case, improved understanding of water clustering on different types of surfaces will be beneficial in studies of various industrial and environmental processes. Supporting information: Figures of experimental water adsorption and desorption isotherms plotted in linear form; figures of experimental nitrogen adsorption and desorption isotherms; a schematic figure of water adsorption and desorption on kaolinite; description of the modelling of water adsorption and desorption on kaolinite. Author contributions A.L and Y.V. designed the study. Y.V. performed experiments. Y.V. and A.A.P. analyzed the experimental data. M.L. performed the molecular dynamics simulations. A.L. and A.N. helped in the conceptual setup of the MD simulations and their interpretation. A.L. performed the adsorption model calculations. A.L., M.L. and A.A.P. wrote the manuscript. A.L. and A.A.P. prepared the figures. All authors discussed the results and commented on the manuscript. Acknowledgements This work was supported by the Research Council of Finland Flagship ACCC (grant no. 337552) and MEDICEN project (grants no. 336557 and 345125). Competing interests The authors declare no competing interests Table of contents graphic