Phase behavior of water-like models in nanoscopic pores of slit shape. Predictions from a density functional theory

We have explored the phase behavior of a set of water-like models in slit pores of nanoscopic dimensions. The interaction between water and pore walls mimics the graphite surface. A version of density functional method is used as theoretical tools. The fluid models are adopted from the work of Clark et al. [Mol. Phys., 2006, 104, 3561]. They reproduce the bulk water vapor-liquid coexistence envelope adequately. Our principal focus is on changes of topology of the phase diagram of confined water and establishing trends of behavior of the crossover temperature between condensation and evaporation on the strength of water-graphite interaction potential. Growth of the water film on the pore walls is illustrated in terms of the density profiles. Theoretical results are discussed in context of computer simulation findings for water models in pores.


Introduction
This manuscript has been prepared as a tribute to Prof. Yuri Kalyuzhnyi, a distinguished scientist in the field of statistical physics and theory of liquids on the occasion of his 70th birthday. Prof. Kalyuzhnyi has made several important contributions to the statistical theory of associating fluids. We appreciate his friendship during last decades. On the scientific side, we have benefited from fruitful discussions with him in the mid of 1990s when we just started the development of the theory of inhomogeneous associating fluids with late Doug Henderson in Mexico [1][2][3]. This brief stage of our close collaboration, during Yu. Kalyuzhnyi's stay at UNAM in Mexico resulted in some novel theoretical insights documented in [4].
Graphitic materials are of much practical importance in several applications that involve solidfluid interfaces (adsorbents, electrodes, solid lubricants). Thus, efficient and precise elucidation of the surface properties of graphite is necessary. Water wettability is one of the most important phenomena characterizing these surfaces. However, from the general perspective, understanding the phase behavior of confined associating fluids and mixtures, which differs significantly from their bulk counterpart, see e.g., [5] for the review, is of undoubtful fundamental and practical importance.
It is difficult to explore the phase behavior of confined fluids by experimental means. Concerning the wetting phenomena of water on graphite or on mica, experiments have been performed using different methodological tools [6][7][8]. Commonly, it is accepted that graphite is hydrophobic. Experimental results usually are given and interpreted in terms of the contact angle values and of the work of adhesion.
Summarizing insight of experimental measurements of the water-graphite system in this aspect is given in table 1 of [7], which provides the list of methods and data for the contact angle.
In general terms, the problem of description of phase behavior of water in different pores using computer simulation techniques is a nontrivial task. Rather comprehensive computer simulation data for the phase behavior of water in cylindrical and slit-like pores of nanoscopic dimensions and with different water-solid interaction strength were reported by Brovchenko et al. [9][10][11][12][13][14] using Gibbs ensemble Monte Carlo methodology for ST2, TIP4P, TIP5P and SPC/E water models. On the other hand, LAMMPS MD method was applied in [15] to describe liquid-vapor phase coexistence of TIP4P2005 water model in graphite and mica slit-like pores.
The most important theoretical findings in this area of research follow from the density functional approach for associating fluids. General aspects of this kind of approaches have been established some time ago [16][17][18]. A summarizing detailed description of the methodology and quite comprehensive insights concerning the results obtained can be found in more recent works [19,20]. For purposes of the present work, we would like to mention that the surface phase transitions in the context of wetting were discussed in [21] for the one-site associating fluid model. The models with two and four-associating sites were studied as well [22,23]. However, the only contribution describing the calculations of the contact angle is [24], up to our best knowledge. Usually, theoreticians resort to the Lennard-Jones potential to describe the non-associating part of the inter-particle interactions. On the other hand, the square-well model has become overwhelmingly most popular with the design of the so-called statistical associating fluid theory (SAFT), see e.g., [25,26] and its application for the modelling of water in [27]. The principal idea behind the modelling of Clark et al. [27] is to reproduce the liquid-vapor (LV) coexistence of water using square-well attraction and site-site chemical association without resorting to electrostatic inter-particle interactions.
In this work, we present a continuation of the project, started in [28,29], focused on the study of the behavior of water in porous media. However, the principal issues considered in the present work have not been explored so far. They are, the surface phase transitions in water-like models confined to wide slit-like pores and their relation to the wetting of graphite by water.

Generalities
We consider a one-component fluid model of associating molecules. Each fluid molecule has four associative sites designated by A, B, C and D, inscribed into a spherical core [30,31]. The set of all the sites is denoted by Γ. The pair intermolecular potential between molecules 1 and 2 depends on the center-to-center distance and orientations, where is the vector connecting site on molecule 1 with site on molecule 2, 12 = |r 12 | is the distance between centers of molecules 1 and 2, is the orientation of the molecule , d is the vector from the molecular center to the site , see also figure 1 of [30]. Each of the off-center attraction sites is located at a distance from the particles' center, = |d | ( = , , , ). In the model in question, only the site-site association AC, BC, AD, and BD is allowed, all association energies are assumed equal. Specifically, the interaction between sites is given as, where as is the depth of the association energy well and is the cut-off range of the associative interaction. The model was used in several previous works.
The non-associative part of the pair potential, ( ), is given as,

33601-2
where hs, ( ) and att, ( ) are the hard-sphere (hs) and attractive (att) pair interaction potential, respectively. The hs term is, hs, where is the hs diameter. The attractive interaction is described by the SW potential, att, where and are the depth and the range of the non-associative attraction potential, respectively.

Comments on the bulk fluid model and its properties
In this work our interest is in four water-like model each with four associating sites designated as W1, W2, W3 and W4 in [27]. In order to obtain a set of parameters to describe density-temperature projection of the phase diagram of liquid water. the following procedure within the SAFT scheme is applied. The free energy of the system is written as a sum of contributions coming from ideal term and residual terms from monomer-monomer non-associative interaction and site-site association, The monomer-monomer term is taken from a Barker-Henderson high-temperature perturbation expansion up to the second order with respect to a hard-sphere reference, whereas the association contribution comes out from the first-order thermodynamic theory of Wertheim [32][33][34][35] in terms of the fraction of molecules not bonded at a site , , All the details concerning these expressions can be found in [27,30,31] and are not given here to avoid unnecessary repetition. In order to obtain five unknown parameters, , , , as and (all defined in the previous subsection), Clark et al. used theoretical predictions for the vapor pressure and saturated liquid density from equations (1)-(3) to reproduce the experimental vapor-liquid equilibrium data for water. The fitting procedure is carried out from the triple point temperature up to the 0.9 fraction of the bulk critical temperature . Due to the degeneracy of the parameters yielding a successful description of the vapor-liquid equlibrium, four optimum sets are proposed, see table 1. They can be referred to as the models with low-dispersion (high hydrogen bonding) values, i.e., W3 and W4, and the models with high-dispersion (low hydrogen bonding) values, i.e., W1 and W2, according to the respective weight of the effects coming from different type of interactions.
Having these parameters available, we employed a simpler version of the theory for the bulk models with the aim to incorporate the necessary expressions into the density functional methodology for inhomogeneous associating fluid. In brief, we used the following treatment [28]. The free energy density (per unit volume), = / is written as, with the following contributions, 33601-3 Table 1. Optimal parameters for water-like model fluids: the diameter of particles, , the depth and range of the square well non-associative potential, and , the cut-off distance of the attractive site-site potential, , and the association energy, as , see [27]. in other words, the hard sphere contribution results from the Carnahan-Starling equation of state and the attractive contribution follows from the mean-field approximation. The notation is as follows, = / , = π 3 /6, Λ is the de Broglie thermal wavelength. Finally, the association term is considered at the level of TPT1 [30,31], with, and where lengthy equation for the site-site bonding volume, , is given for example by equation (20) of [28]. The expressions for the chemical potential, , and pressure, , involved to study vapor-liquid equilibrium are omitted for the sake of brevity and to avoid unnecessary repetition, they are given by equation (22) and (23) of [28]. The dimensionless input parameters for all the models in questions are given in table 2. In the same table, for the sake of convenience of the reader, we provide the values for the critical temperature that results from our version of the theoretical procedure for each of the bulk models. Concerning dimensionless units used below, we introduce them in the following common form, * = 3 , * = / , * = / . Table 2. Dimensionless parameters for the W1, W2, W3 and W4 four-site water-like models: the association energy, * as = as / , the site-site bonding volume, * = / 3 , the cut-off distance of the attractive site-site potential * = / ( / = 0.25). Critical temperature of the models, * , is given in the last column. Two projections, * - * and - * , of the phase diagrams for the bulk models in question are shown in figure 1. Both are in reduced units, i.e., rescaled to the respective values at a critical point. Somewhat unexpectedly, one can observe that the models behave according to the law of corresponding states. Actually, it means that in spite of quite different values of parameters for the square-well attraction and for the association energy, the balance between these two effects is preserved for all four models in question.

Model
Two defficiencies of the applied equation of state should be mentioned. One of them, mentioned already by Clark et al. [27], is that water anomalies are not captured by the present modelling. On the other hand, our preliminary calculations have shown that the critical compressibility factor of this set of models is not satisfactory. These issues should be analyzed in more detail having in mind the available analyses, see [36,37]. The dimensionless critical temperature for each of the models from table 2 can be converted into real units to yield: = 679.5 K for W1, = 659 K for W2; = 585.2 K for W3 and finally = 508.6 K for W4 model. Thus, the performance of the W1 and W2 models in this aspect is much better compared to other two models (the experimental result is at 647 K). The critical temperature for the most popular water models is reported in table 2 of [38]. It ranges from 640 K for TIP4P/2005 and 638.6 K for SPC/E models down to 521 K for TIP5P water model. In addition to the fraction of nonbonded species at a site, it is possible to evaluate the fraction of particles in different bonding states [30,31]. Previous calculations of these fractions for the model with four associating sites referred to as associating hard spheres [39,40]. Here, we would like to show the fractions of differently bonded species along the bulk coexistence envelope for all four models in question. It appears that in spite of quite different balance between non-associative attraction and site-site bonding, the behavior of ( * ) (on reduced temperature to make comparison reasonable) is quite similar. Minor differences are observed between the performance of W1 and W2 models, whereas the W3 and W4 models lead to practically indistinguishable results. The vapor phase  is predominantly composed of nonbonded species. In liquid phase, however, the maximum magnitude of fractions with one, two, three and four bonds is approximately the same. At a high temperature, the dominant fractions are 1 and 2 , whereas at a low temperature (of the order of possible triple point), the fluid structure is dominated by 3 and 4 describing well defined hydrogen bonding network. From theoretical perspective, it seems worth to study the application of a more sophisticated procedure to describe the non-associative attraction, e.g., the first-order mean spherical approximation (a limited set of results are reported in [28]), and to test how sensitive is the association contribution to the description of the reference hard-sphere fluid. These issues require a separate investigation.

Fluid-wall interaction
The pore walls are located at = − /2 and = /2. The external potential, ( ), exerted on a fluid particle inside the pore by walls is, The function ( ) is given by the Steele's 10-4-3 gas-solid potential [41,42], where sf , sf are the cross energy and the size parameters describing solid-fluid interaction. The characteristics of graphite are as common, Δ = 0.335 nm, = 114 nm −3 , ss / = 28 K, ss = 0.34 nm. One frequently used possibility is to evaluate the cross interaction parameters sf and sf by assuming the Lorentz-Berthelot (LB) combination rules, sf = √ ss , and = ( ss + )/2. However, one can consider these two parameters as free as well. In what follows, we use the abbreviation for the energy of interaction of water-like species with the pore walls denoted as * gs , * gs = 2π sf 2 sf Δ/ . The pore width is given in dimensionless units as well, * = / .

Theory
The system is studied using a version of the density functional theory (DF), described already in detail in [29,43,44]. To avoid unnecessary repetition, all this stuff is given as appendix. Rather, we proceed to the description of the results and novel observations.

Results and discussion: phase diagrams, adsorption isotherm, density profiles
The entire set of water-like models of this study is considered under confinement in a wide slit-like pore, * = 21, to describe surface phase transitions and respective phase diagrams. It is worth mentioning that the computer simulation results of Brovchenko et al. [9][10][11][12][13][14] are obtained for narrower pores, in the case of the slit shape pores the maximum width is of the order of * ≈ 9.5. In order to cover an ample set of effects due to changes of the strength of water-pore walls interaction, these authors assumed different values for the parameter gs (ranging from −0.39 kcal/mol to −4.62 kcal/mol) to mimic hydrophobic paraffin-like substrate, moderately hydrophilic, carbon-like surface, and strongly hydrophilic, silica-like surface. On the other hand, Srivastava et al. [15] refer to the pores of the maximum width, * = 19, to describe hydrophilic mica pore and less hydrophilic pore with graphite walls. It is worth noting that the unique 9-3 water-substrate interaction potential was assumed in the works from I. Brovchenko laboratory, whereas in [15] the potentials for graphite and mica pores differ from each other.
To keep track with these computer simulation results, in the present work we employ the following nomenclature of pores: if the value of * gs parameter is small, so that capillary evaporation is observed in the entire range of temperature, the pore is termed as hydrophobic. By contrast, if capillary condensation  is observed in the entire temperature range for big * gs values, the pore is termed as hydrophilic. The intermediate situation, with condensation at high temperatures and evaporation at low temperatures, is termed as moderately hydrophilic pore. Let us proceed now to the results for water-like models under different confinement.

Hydrophobic pores
To begin with, we consider the models W1 and W4 in the hydrophobic environment. The first of them, W1 model, is characterized by a rather wide square-well for nonassociative interparticle attraction, whereas the W4 model has the most narrow attractive well but the most strong association interaction, cf. table 1. Different projections of the phase diagrams are given in figure 3.
Principal theoretical observation is that for both models confined in hydrophobic pore, the liquiddensity branch of the coexistence envelope pronouncedly departs from its bulk counterpart. Magnitude of this effect depends on the value of * gs (figure 3a). The vapor branch of the coexistence envelope is much less affected by confinement. The critical temperature decreases due to confinement and the critical density for both models in the pore is slightly smaller than the bulk critical density * . These DF predictions are in agreement with the trends observed in computer simulation studies, cf. figure 1 of [14] and figure 2 of [12]. In the case of a wide graphite-like pore, the critical density is almost the same as the bulk value, cf. figure 2a in [15] but depletion of the critical temperature is well marked. The shift of the critical density is correlated with the depletion of the liquid density branch for confined model w.r.t. its bulk counterpart.
Changes of water structure upon capillary evaporation in hydrophobic pore are illustrated in figure 4. The principal effect is that the fluid density profile is strongly depleted close to the pore walls. The  magnitude of depletion depends on * gs and on temperature. If the density profiles for two models are plotted at the same reduced temperature, * = * / * , cf. panels a and b of figure 4 at the same * gs , the distribution of molecules in the pore is very similar. In qualitative terms, the shape of the density profiles coming from density functional theory and from simulation data cited just above is quite similar.

Moderately hydrophilic and hydrophilic pore
The phase behavior of water models in hydrophilic pores is much richer in comparison to hydrophobic pores. Our first comment refers to the presentation of some plots below. In several previous works from this laboratory, e.g., [18,22,29] we have used the 0 / - * plots for the phase diagrams [ 0 ( * ) is the chemical potential of vapor-liquid bulk coexistence at temperature * ]. On the other hand, the ( − 0 ) versus * projection is formally equivalent but it is more convenient to distinguish between condensation and evaporation in the pore and most importantly to establish if the system approaches the wetting borderline or departs from it. Moreover, this type of presentation was used in the classical work [45] to perform classification of the surface phase transitions based on the strength of attraction between fluid and substrate.
The curves in figure 5 illustrate the phase diagrams of W1 model for different values of the watersubstrate attraction in terms of * gs . The case corresponding to * gs = 8.311 that follows from the LB combination rules is marked as CR in the figure.
In contrast to the discussion concerned with the behavior of water in hydrophobic pore, the phase diagrams in figure 5a exhibit shrinking of the coexistence envelope under confinement that involves vapor and liquid branches w.r.t. bulk coexistence. Depression of the critical temperature is observed. However, the critical density for a confined model now is higher than in the bulk. The vapor branch can be a single-valued function ( * gs = 6) or split due to the appearance of additional phase transition ( * gs = 9.5). This is illustrated even better on the chemical potential -temperature projection, figure 5b. For moderately hydrophobic pore at relatively low values of * gs , the phase diagram is a single branch that takes negative values along − 0 axis at high temperatures, i.e., the fluid vapor converts into liquid in the entire pore prior bulk transition (capillary condensation). At low temperatures, the pore filling occurs at positive values of − 0 , i.e., capillary evaporation is observed. The crossing point for each * gs value under such circumstances describes transition between two regimes. Thus, the coexistence envelope plotted in figure 5a for the case * gs = 6, describes condensation in the pore for * > 1.7 and capillary evaporation at lower temperatures.
If the value of * gs is higher, one observes two scenarios. The phase diagram may be a single branch, e.g., at * gs = 8.311, or it consists of more branches, e.g., * gs = 9.5 (figure 5b). In the first case, metastable layering is observed whereas in the second case stable part of the diagram at low temperatures describes the formation of the fluid layer on the pore walls. This part joins the principal part of the coexistence envelope at a triple point. At still lower temperatures, the formation of such a layer is metastable w.r.t. the entire pore filling or capillary condensation. Moreover, the inclination of the − 0 - * projection of the phase diagram at low temperatures is different (figure 5b). It either tends to the = 0 line or departs from the bulk coexistence in the entire temperature interval studied.
The results from figures 5b do not permit to get deeper insights into the entire vapor side behavior of coexistence and in particular to establish what kind of phase is formed upon layering transition. To do that, one should return to - * projection, figure 5a. First interesting observation concerns reasonably high temperatures. At a low value of * gs , e.g., * gs = 6 (black lines), vapor phase is quite "dilute", whereas at * gs = 8.311, the "dilute" phase is approximately twice denser than prior to transition to the entire pore filling. Anyway, the water film continuously grows prior to condensation transition. The transition can be interpreted as thin film (vapor) -filled pore ( * gs = 6) or thicker film (still vapor) -filled pore transition ( * gs = 8.311). However, at even stronger fluid-substrate attraction, the growth of the fluid film on the pore walls from the "underlying vapor" occurs via discontinuous formation of a separate layer phase that subsequently converts discontinuously into liquid upon condensation in the pore. At even lower temperatures, vapor converts into liquid discontinuously but without an intermediate phase. Thus, there exists a certain interval of * gs and the corresponding temperature interval in which a water film grows either continuously or discontinuously prior to filling the entire pore. Such changes upon increasing the 33601-9 chemical potential are accompanied by changes of the fraction of the average number of non-bonded at a site molecules. This kind of projections of the phase diagram are shown in figure 5c. On the vapor side, the fraction of non-bonded species is very high, whereas the liquid phase structure is characterized by a much higher fraction of species bonded between themselves due to site-site association.
It is worth to comment our observations in view of computer simulation results. Namely, the situation when the critical density for the confined fluid becomes higher than for the bulk due to an increasing hydrophilicity was observed and discussed in [15], cf. figures 2a and 2b of [15]. Nevertheless, a layering transition was not observed with their model of fluid-substrate interaction. Computer simulation evidence of changes of the vapor branch of the coexistence envelope in slit-like hydrophilic pores was documented in [12] (figure 5) and in [13] (figure 10). These results refer to the ST2 water model. The authors consider that the formation of an additional phase is attributed to prewetting. We do not share this opinion completely because a complementary exploration of one-wall phenomena such as wetting for the models in question is necessary to make definite conclusions.
Deeper and complementary insights for the present discussion should come out from the exploration of adsorption isotherms and fluid structure in the pore. The phase diagrams in figure 5 were constructed from a large set of adsorption isotherms. For purposes of analysis, in figure 6 we show a comparison of adsorption isotherms obtained at a fixed temperature * = 2.4, and at two values for the strength of water-substrate interaction, cf. figure 5a. At a higher value of * gs , * gs = 9.5, the water film on the walls grows faster and the average density in a wide pore before transition is approximately twice higher than for the case of * gs = 6. More interesting is the evolution of the density profiles of molecules shown at certain characteristic points of the isotherm with * gs = 9.5. First, a layer adjacent to the pore walls is formed (at a point 1). The height of the first maximum increases if one "moves" to point 2. However, the growth of the first layer is accompanied by a simultaneous formation of the second layer, due to inter-particle attraction (recall that the square well attraction is rather wide for this model, cf. table 1) and presumably due to inter-particle bonding. At point 3, just before pore filling, both the first and the second layer are characterized by quite high maxima of the density profile, so that they are almost entirely filled. There is a signature of the third layer as well. Upon condensation, the already well formed first and second layers do not suffer big changes, in contrast to the third layer. In the entire space, farther from the pore walls, the vapor drastically converts into liquid.
In order to elucidate more details of the changes of the structure, we compared the adsorption isotherms of the same system with * gs = 9.5, but at two temperatures, figure 7a. At a lower temperature, * = 2.15, the first-order layering transition occurs (points 1 and 2). At a slightly higher value of the chemical potential, the pore filling (capillary condensation) takes place. The metastable portion of the adsorption isotherm is shown in figure 7a to illustrate that the growth of the film is smooth, without layer by layer filling. The density profiles in figure 7b describe how the water structure changes along the isotherm. The first layer, adjacent to the wall, is quite dense even before layering transition. After layering transition, the first layer almost reaches its maximum capacity. The principal change is the formation of the second layer. There is no signature of the third layer. This kind of "bilayer" structure converts into a liquid structure by completion of the second layer together with the filling of the entire pore.
Proceed now to the description of the results for other water-like models. According to the parameters from table 1, one would expect that the phase behavior of the W2 model exhibits a similarity to the one for W1 model. This appears to be true. However, some details are new, in comparison to the findings for the W1. The phase behavior of the W2 model in a wide pore at different strengths of watersubstrate interaction is shown in figure 8. In figure 8a, we observe that if * gs is chosen from the LB combination rules, stable layering transition is observed at certain temperatures, in contrast to W1 with * gs chosen according to combination rules. This transition describes the formation of "bilayer" structure. Moreover, the layering transition envelope joins the principal coexistence at the triple point. Under a further increase of * gs , the phase diagram exhibits interesting changes. The layering 1 + 2 transition in the upper part, decouples into two true, layer by layer, transitions. Both of them join at the triple point. The critical temperature of the first layering is lower than of the second one. At lower temperature, when the effects of inter-particle and of particle-wall attraction become stronger, water again forms a film with two layers simultaneously. Therefore, decoupling of the upper part should be attributed to weakening of inter-particle, presumably non-associative, attraction. These events are observed in the system when the coexistence line departs quite far from the bulk, see blue lines in figure 8b.
Insights into the changes of water structure close to the pore walls are illustrated in figure 9. The sequence of changes of the profiles exhibits a similarity to the trends discussed for figure 7b. However, it is worth mentioning that the first layer, already after the second layering transition reaches its maximum capacity. The second layer is also well filled at these conditions and there is no signature of the third layer prior to filling of the entire pore, figure 9a.
Additional view into the structure is provided in figure 9b, where we plot the density profile for particles that are not bonded at a site. The fraction ( ) follows the trends of the behavior of density  in figure 9a. Higher local density results in lower local values for the fraction of non-bonded species or equivalently causes a stronger bonding tendency. In this plot we observe that after the second layering transition, the bonding state within the first layer very closely resembles the bonding state of the liquid phase filling the pore. The second layer is characterized by a high degree of bonding as well. How good are these values in particular states would require performing quite elaborate computer simulations. More pronounced peculiarities of the phase behavior are expected while considering the W3 and W4 models in comparison to the W1 and W2 ones. The principal reason is that the former two models involve a much narrower square attractive well, cf. table 1. Immediate consequences can be seen from the phase diagrams of the W3 model in figure 10. Even for a quite weak fluid-substrate attraction, e.g., * gs = 5, the pore wall is capable of pulling the water molecules and forming a "bilayer" structure. However, this piece of the coexistence is separated from the principal part. For lower temperatures, water vapor converts into liquid without any peculiarities. This kind of separation results in a phase diagram consisting of three separate branches as shown in figure 10b. The low temperature branch is just the vapor-liquid coexistence in the pore, i.e., the capillary condensation. Two branches at higher temperature describe a consecutive layering (1 + 2 branch) and subsequent condensation that starts from a thin bilayer film. This thin film is even wider, if one focuses into even higher temperatures. This behavior differs from the trends observed for W1 and W2 models discussed above.
If the strength of water-substrate interaction is chosen from the combination rules, the phase diagram becomes even more complex. Namely, the pore wall is capable of decoupling the 1 + 2 "bilayer" film into two consecutive layering transitions in the upper part of this fragment of the diagram with the triple point between them. The critical temperature of the first layering transition is higher than of the second one, in contrast to what we have seen just above, describing the phase behavior of the W2 model. Returning to figure 10a, we observe that at lower temperatures, the formation of the "bilayer" structure via a single transition is followed by one more transition describing the formation of the third layer. Next, this multilayer structure converts into liquid upon increasing the chemical potential. However, the fragments of the phase digram describing the formation of layers are not connected to the principal envelope, i.e., there are no triple points. A better illustration of this behavior is shown in figure 10b by red lines. The only triple point is seen between the first and second layering transitions.
More complexity is expected for the W4 model since its parameters correspond to the most narrow square well compared to other models. The results of our calculations illustrate this in two panels of figure 11. The phase diagram evolution with increasing * gs from 2 to 3 keeps track of our discussion while moving from W2 to W3 model. In both cases, a "bilayer" structure is formed at higher temperatures, whereas at lower temperatures, capillary condensation occurs.
Next, the phase diagram corresponding to * gs = 4 exhibits a similarity to the phase behavior of the W3 model. However, we observe an enhanced difference between the critical temperatures of the first and second layering transitions that join at the triple point. These structures then transform, because of the additional transition due to the formation of the third layer. Afterwards, capillary condensation occurs upon increasing the chemical potential as expected. The final result in this figure concerns the * gs value arising from the combination rules. None of the triple points is seen. All three layering transitions are separated between themselves and separated from the principal condensation envelope.
Additional insight into how the water film grows on the pore walls within these models with a narrower square well (W3 and W4) comes out from figure 12. We compare two examples of the adsorption isotherms at the same reduced temperature and in both cases consider the fluid-wall interaction strength value from the combination rules. It can be seen that the parameter (square well width) is of much importance. If the attractive square-well is narrower (at a fixed, sufficiently strong association energy), then the tendency for the formation of several layers on the pore walls is stronger. The metastable parts of the isotherms consist of consecutive layering transitions. These metastable parts can be compared with their entirely smooth "counterpart" for the W1 model (with a much wider square well attraction) in figure 7a. The phase diagrams obtained and discussed above in detail permit to search for * gs values satisfying the = 0 condition at various dimensionless temperatures, * . This search is partially illustrated in figures 5 and 8. Next, with * gs and * available, one can obtain the corresponding sf values. The crossover temperature as a function of water-graphite interaction energy is shown in figure 13.
Two observations can be deduced from these curves. Similarity of the behavior of W1 and W2 models is due to a rather wide square-well attraction of these models. If the balance between the non-associative attraction and associative site-site bonding is different, like in W3 and W4 models, a shift of the curves to lower values of sf is observed. However, the inclination of the curves in the entire set is similar. If sf decreases, the interval of temperature, where one observes a capillary condensation, shrinks or, in other words, an interval of temperature describing evaporation extends. In qualitative terms, this behavior can correspond to changes of the wetting temperature. One can attempt to evaluate the crossover temperature in pores a few times wider than we study at present, in order to approach the description of one-wall phenomena. However, a strict alternative is to set up the problem with a single wall and permit the growth of the "macroscopic" film with free terminating interface.

Summary and conclusions
The water-like models of this study were developed by Clark et al. [27] and describe vapor-liquid coexistence of the bulk water pretty well. The inter-particle interaction potential consists of spherically symmetric hard core, square-well attraction and associative site-site attractive interaction that mimics hydrogen bonding effects. However, we have explored the phase behavior of these models confined to a quite wide slit-like pore. To do that, a version of classical density funcional theory is used. The principal issue we dealt with is in changes of topology of the phase diagrams upon changing the strength of water-substrate interaction. We considered a hydrophobic pore as well as moderately hydrophilic and hydrophilic environment for water species. A similar type of investigation was performed in [9][10][11][12][13][14][15] using computer simulations techniques. In spite of these efforts, still the problem is not comprehensively studied. A theoretical approach of the present work has an advantage over simulations because it intrinsically involves the chemical potential variable indispensable in the study of phase equilibria.
The water-substrate interaction has the form of 10-4-3 potential formulated by Steele [41] and commonly used in adsorption studies on a graphite surface. The combination rules were used while deriving values for energy of interaction between water molecules and graphite. On the other hand, we assumed various values of the energy to explore in detail the phase behavior of water under confinement. Our findings qualitatively agree with computer simulation results. In hydrophobic pores, interpretation of the behavior of water is much easier than in hydrophilic pores. In the latter type of systems, the phase behavior is quite complex and difficult to access by computer simulation methods. For moderately hydrophilic and hydrophilic setup, we established conditions at which additional phases appear at the vapor side of the coexistence envelope. This is done for all four water-like models in question. Stable and metastable parts of these rather low density phases were evaluated. Triple points under different circumstances were obtained. Continuous and discontinuous growth of the water film on substrates with different affinity to water was described. In certain situations, we observed the formation of a bilayer-type water film and the possibility of its decomposition into layering transitions. Exploration of the phase  behavior is complemented by the description of adsorption isotherms and of the density distribution of molecules adjacent to the pore walls. Our results may be combined with the development of studies of adsorption of water in various pores, e.g., [23,46,47].
Concerning the missing elements, we would like to mention that the problem of wetting of this class of models for an associating fluid on graphite requires a separate investigation. This should permit to evaluate the wetting temperatures and the contact angles. Important and interesting extensions should involve exploration of various mixtures in all aspects of the present work, see e.g., [48] as a possible starting point.
where Φ as is the associative free energy density [29,43,44] that results from a generalization of the firstorder thermodynamic perturbation theory (TPT) of Wertheim to nonuniform systems. Similarly to the case of the hard-sphere contribution,˜h s , the evaluation of the functional Φ as (r) requires the knowledge of the weighted density profiles (r) and n (r) [17]. We have [29,43,44] Φ as ( ) = 4 0 (r) (r) ∑︁ =1 ln (r) − 1 2 ( (r) − 1) , (A. 6) where, in contrast to the bulk system, the fraction of non-bonded molecules at a given site , depends now on the position vector r. All the details and the definition of the quantity (r) are given in previous works [29,43,44]. We only note that the evaluation of the functional Φ as (r) requires the knowledge of the contact value of the pair distribution function. The latter is computed from a generalizations of the Carnahan and Starling equation of state to nonuniform systems [17]. The density profile results from the minimization of the Grand Canonical potential, cf. [20], where is the bulk reference system, i.e., the bulk, uniform system at the same temperature and the chemical potential as the system under study. The quantity is the configurational chemical potential of the bulk system = − id . The average density of the fluid in pore, , is calculated from the density profile,