505Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
DOI: 10.17344/acsi.2021.7048
Feature article
Analytical Modeling of Thermodynamic and Transport 
Anomalies of Water
Tomaz Urbic*
University of Ljubljana, Faculty of Chemistry and Chemical Technology, Chair of Physical Chemistry,  
Večna pot 113, SI-1000 Ljubljana, Slovenia
* Corresponding author: E-mail: tomaz.urbic@fkkt.uni-lj.si
Received: 07-12-2021
Abstract
The structures and properties of biomolecules like proteins, nucleic acids, and membranes depend on water. Water is also 
very important in industry. Overall, water is an unusual substance with more than 70 anomalous properties. The under-
standing of water is advancing significantly due to the theoretical and computational modeling. There are different kinds 
of models, models with fine-scale properties and increasing structural detail with increasing computational expense, 
and simple models, which focus on global properties of water like thermodynamics, phase diagram and are less com-
putationally expensive. Simplified models give a better understanding of water in ways that complement more complex 
models. Here, we review analytical modelling of properties of water on different levels, the two- and three-dimensional 
Mercedes– Benz (MB) models of water and experimental water. 
Keywords: Water, statistical model, anomalies, transport properties, analytical model
1. Introduction
The Earth is a watery place by water being the most 
important fluid in nature for life and for humans in the 
industry.1–4 About 71 percent of the Earth’s surface is wa-
ter-covered, and the oceans hold about 96.5 percent of all 
Earth’s water. Water also exists in the air as water vapor, in 
rivers and lakes, in icecaps and glaciers, in the ground as 
soil moisture etc. Water controls the planets geochemical 
cycles; is a dominant driver of biomolecules, drug interac-
tions, and biological actions; and central to green chemis-
try and many industrial processes.5,6 Water is essential for 
our bodies. Every major system in our body depends on 
water to function since approximately two thirds of your 
body is water. Due to all these facts the structure and ther-
modynamics of water and aqueous solutions is of great 
importance for all sciences, especially chemistry and biol-
ogy. Water exhibits many anomalous properties that affect 
life at a larger scale. Many animals benefit from the large 
latent heat of water to cool them down through sweating. 
The large heat capacity of water prevents local temperature 
fluctuations, facilitating thermal regulation of organisms. 
The density anomaly and lower ice (hexagonal ice Ih) den-
sity have a huge effect on surviving of organisms in frozen 
seas and lakes. Water is an almost universal solvent.7 Near-
ly all known chemical substances will dissolve in water at 
least to a small extent. In comparison to other liquids, it 
has the most puzzling behavior.7,8 It is said that water is an 
anomalous liquid.  Anomalous liquids are liquids that ex-
hibit unexpected behavior upon variations of the thermo-
dynamic conditions in comparison to normal (argon-like) 
liquids. Water is the classic example of such anomalous 
liquids. Water’s density maximum at 4 °C, the lower densi-
ty of the solid phase compared to the liquid phase, high 
and nearly constant heat capacity in the liquid phase, neg-
ative expansion coefficient below the temperature of the 
density maximum, as well as high surface tension and vis-
cosity are the most known examples of anomalous proper-
ties. If we continue, we have an anomalous increase in the 
compressibility and specific heat by cooling, unusually 
high boiling, freezing and critical points. The reason for 
water’s complexity is due to its strong orientation-depend-
ent hydrogen bonding and strong intermolecular associa-
tions.
An understanding of the hydrogen bonding is there-
fore crucial to understand the behavior and properties of 
water and aqueous solutions. Yet, despite extensive theo-
ry and simulations, the fact how water’s properties come 
from its molecular structure remains poorly understood. 
Many models of varying complexity have been developed 
and analyzed to model water’s extraordinary properties, 
506 Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
for review see literature.7–19 The rigid models are consid-
ered the simplest water models and rely on non-bonded 
interactions. The electrostatic interaction is modeled using 
Coulomb’s law, and the dispersion and repulsion forces us-
ing the Lennard-Jones potential. Examples of such models 
are SPC (simple point-charge)20, TIP3P (transferable in-
termolecular potential with 3 points)20 and TIP4P (trans-
ferable intermolecular potential with 4 points)21 etc. In po-
larizable models we consider many-body energies which 
can be effectively accounted for by a single term represent-
ing classical many-body polarization. Several polarizable 
water models, with different degrees of sophistication, 
have been developed and used in molecular dynamics 
and Monte Carlo simulations of aqueous systems.22 A key 
goal of the liquid-state statistical thermodynamics is to 
develop a quantitative theory for water and aqueous solu-
tions. Theory and simulations have not yet been able to 
explain how water’s molecular structure leads to its den-
sity, compressibility, expansion coefficient and heat capac-
ity as functions of temperature and pressure, including 
its well-known anomalies. The properties of water can be 
determined with quantum-mechanical calculations.22,23 
These methods offer the highest degree of exactness, but 
a high computational cost of these approaches limits their 
use to small water systems, even though these insights 
allow the development and fine-tuning of simplified wa-
ter models.24–26 There have been two main approaches to 
modeling liquids. One approach is to perform computer 
simulations of atomically detailed models. Another way 
captures many properties of water and aqueous solutions 
by simpler models. 
One of the simplest models for water is the so-
called Mercedes-Benz (MB) model,27 which is a 2-di-
mensional model that was originally proposed by Ben- 
Naim in 1971.28,29 Each MB water particle is modeled 
as a disk that interacts with other particles through: (1) 
a Lennard-Jones (LJ) interaction and (2) an orienta-
tion-dependent hydrogen bonding interaction through 
three radial arms arranged as in the MB logo. Interest in 
simplified models is due to insights that are not obtaina-
ble from all-atom computer simulations. Simpler models 
are more flexible in providing insights and illuminating 
concepts, and they do not require big computer resourc-
es. The analytical models can also provide functional re-
lationships for engineering applications and lead to im-
proved models of greater computational efficiency. For 
the MB model, the NPT Monte Carlo simulations have 
shown that it predicts qualitatively the density anomaly, 
the minimum in the isothermal compressibility as a func-
tion of temperature, the large heat capacity, as well as the 
experimental trends for the thermodynamic properties 
of solvation of nonpolar solutes27 and cold denaturation 
of proteins.30 The MB model was also extensively stud-
ied with analytical methods like integral equation and 
thermodynamic perturbation theory31–36 and statistical 
mechanic modeling37–39. Recently also phase diagram of 
liquid part and percolation curve of the model was calcu-
lated and reported.40 The MB model has also been used 
to study systems with water molecules confined in par-
tially quenched disordered matrix41–43 and within small 
geometric spaces.44,45 Nonequilibrium Monte Carlo and 
molecular dynamics simulations were used to study the 
effect of translational and rotational degrees of freedom 
on the structural and thermodynamic properties of this 
MB model.46–48 By holding one of the temperatures 
constant and varying the other one, the effect of faster 
motion in the corresponding degrees of freedom on the 
properties of the simple water model was investigated. 
The situation where the rotational temperature exceed-
ed the translational one is mimicking the effects of mi-
crowaves on the water model. A decrease of rotational 
temperature leads to the higher structural order while an 
increase causes the structure to be more Lennard-Jones 
fluid like. The 2D MB model was also extended to 3D by 
Dias et al.49 and Bizjak et al.50,51 
Even though computer simulations play an impor-
tant role in understanding the properties of liquids, they 
can be quite time consuming, even for simple two-dimen-
sional 2D models. Due to this it is equally important to 
develop simplified, more analytical approaches. One such 
model is a statistical mechanical model, developed by 
Urbic and Dill33. The model is directly descendant from 
a treatment of Truskett and Dill, who developed a nearly 
analytical version of the 2D MB model.52,53 In the model, 
each water particle interacts with its neighboring parti-
cle through a van der Waals interaction and an orienta-
tion-dependent interaction that models hydrogen bonds. 
Recently this theory was extended to 3D MB model38 and 
later parametrized to describe properties of experimental 
water.54 
In this paper, we made review of analytical modeling 
for MB model of water, its properties in bulk which are 
starting point to develop the theory for solvation of polar 
and nonpolar solutes, important for example in self-as-
sociation of surface-active compounds such as ionic liq-
uids,55 protein folding, etc. The outline of the paper is as 
follows. We present the 2D and 3D MB model in Sec. 2, 
and the details of the statistical mechanical methods are 
done in Sec. 3. In Sec. 4 we show and discuss the results 
and summarize everything in Sec. 5.
2. The Model
2. 1. 2D MB Model
In 2D, the water particles are modelled as a two-di-
mensional disk with three bonding arms separated by an 
angle of 120°, which is fixed as in Mercedes-Benz logo (See 
Figure 1).27 These arms mimic formation of hydrogen 
bonds. The interaction potential between particles i and j 
is a sum of a Lennard-Jones (LJ) and a hydrogen-bonding 
(HB) term
507Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
 (1)
Where rij is the distance between centers of parti-
cles i and j. Xu i, X
u 
j are the vectors representing the coordi-
nates and the orientation of the particles i and j. The Len-
nard-Jones part has a standard form
 (2)
sLJ and eLJ are the depth and the contact value of the LJ 
potential. The hydrogen bonding part is the sum of inter-
actions between all pairs of the arms of different molecules
 (3)
and is described by Gaussian function in distance and 
both angles
 (4)
Here, eHB = –1 is a HB energy parameter and rHB = i 
is a characteristic length of HB, u ij is the unit vector along 
ru ij and iuk is the unit vector of the kth arm of the particle I, 
and qi is the unit vector of the ith arm of the particle j. qi, qj 
are the orientations of the particle with respect to x axes. 
G(x) is the unnormalized Gaussian function
 (5)
The strongest hydrogen bond occurs when an arm 
of one particle is co-linear with the arm of another par-
ticle and the two arms point in opposing directions. The 
LJ well-depth eLJ is 0.1 times the HB interaction energy 
|eHB| and the Lennard-Jones contact parameter sLJ is 0.7 
rHB. The width of the Gaussian function for distances and 
angles (s = 0.085 rHB ) is small enough that a direct hydro-
gen bond is more favorable than a bifurcated one.
2. 2. 3D MB Model
In 3D, each water molecule is represented as a Len-
nard-Jones sphere (LJ) with four arms oriented tetrahe-
drally.50 The angle between neighboring arms in a mole-
cule is 109.47° (see Figure 2). Like in 2D, in 3D the inter-
Figure 1: The MB particles in 2D.
Figure 2: The MB particles in 3D.
508 Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
action potential between two water molecules is a sum of 
the Lennard-Jones potential and the hydrogen bond term
 (6)
The Lennard-Jones part of the potential is the same 
as in 2D. The hydrogen bonding part of the interaction po-
tential is
 (7)
Where Ωu i. Ω
u 
j are the orientational vectors of both 
particles and UklHB (rij, Ω
u 
i, Ω
u 
j) describes the interaction 
between two HB arms of different molecules
 (8)
Like in 2D, the strongest hydrogen bond occurs 
when an arm of one particle is colinear with the arm of 
another particle pointed towards each other. The mod-
el does not make a distinction between hydrogen bond 
donors and acceptors. Apart from the dimensionality, we 
want to keep the 3D MB model as similar as possible to the 
original 2D MB model. Hence, the parameters of our 3D 
model are the same as used in the 2D MB model calcula-
tions, except for the depth of the Lennard-Jones potential 
well eLJ. This change was made to maintain the same ratio 
between strength of the Lennard-Jones interaction and hy-
drogen bond interaction due to the different geometries; 
eLJ=1/35 eHB. These model parameters were not chosen or 
optimized to compare with experiments and can undoubt-
edly be improved for those purposes.
3. The Statistical Mechanics Theory
3.1. 2D MB Model
In the theory, the system consists of N water mol-
ecules.37 To keep track of the state of interaction of each 
possible hydrogen bonding arm of each water molecule we 
are using an underlying ice lattice as a bookkeeping tool. 
For the 2D water model, the underlying lattice is hexago-
nal (See Figure 3). We focus on a single water molecule in 
the hexagon and the relationship of that water to its clock-
wise neighbor. Figure 4 shows the three possible relation-
ships: the test water can either form a hydrogen bond, a 
van der Waals contact, or no interaction at all. We compute 
the isothermal-isobaric statistical weights, ΔHB of the hy-
drogen-bonded molecules, ΔLJ of the van der Waals con-
tacts, and Δ0 of the unbonded population as functions of 
temperature, pressure, and interaction energies.
 The hydrogen-bonded state. If the test water mole-
cule points one of its three hydrogen bonding arms at an 
angle θ to within π/3 of the center of its clockwise neighbor 
water, it forms a hydrogen bond. The energy of interaction 
of the test water is
 (9)
k is the angular spring constant that describes the weak-
ening of the hydrogen bond as it becomes increasingly 
off-angle, and eHB and eLJ are the potential energy parame-
ters. We regard this type of hydrogen bond as weak insofar 
as it is not cooperative with neighboring hydrogen bonds. 
We consider a more cooperative type of hydrogen bonding 
below. To compute the isothermal-isobaric partition func-
Figure 3: The lattice of the model showing both the hexagon of the icelike structure and illustrating a pair interaction used for bookkeeping to avoid 
triple counting.
509Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
tion, ΔHB, of this HB state, we integrate this Boltzmann 
factor over all the allowable angles and over all the allow-
able separations x and y of the test molecule relative to its 
clockwise neighbor,
 (10)
Where b = 1/kBT is inverse temperature, c(T) is the 
2D version of the kinetic energy contribution to the parti-
tion function and vHB is volume per molecule in this state. 
òò dxdy represents the volume over which the second 
molecule has translational freedom to form a hydrogen 
bond with the first water and is equal to the effective vol-
ume . The volume vHB of the hydrogen-bonded 
state is determined in the following way. First, we estimate 
an upper bound on the volume, from a simple geometric 
calculation. For the perfect hexagon crystal, representing 
low-pressure ice, the volume of the solid if the center of the 
hexagon is unoccupied is
 (11)
Second, since liquid water is denser than ice, we es-
timate a lower bound on the volume using high-pressure 
ice, where another MB water occupies the center of each 
hexagonal cage52,53
 (12)
Since the density of liquid water must lie between 
these limits, we estimate its volume as
 (13)
Where xv = 1.01 is chosen empirically by fitting the 
density dependence vs. temperature. Using these defi-
nitions and performing the integration in Equation (10) 
gives
 (14)
The van der Waals (vdW) state. Here, the test water 
molecule forms only a van der Waals contact with its clock-
wise neighboring water. The water molecule has an energy
 (15)
The isothermal-isobaric partition function, ΔLJ of 
this state is given by integrating over angles and positions 
of the test particle relative to its clockwise neighbor as in 
case of the HB state
 (16)
The integral òò dxfy represents the translation vol-
ume over which the second molecule forms a van der 
Waals contact with the first water and is equal to the effec-
tive volume veffLJ = 0.104. The volume occupied by mole-
cule in this state, vLJ, is volume of packed LJ disks
 (17)
Integrating the partition function gives 
 (18)
Figure 4: The three model states: (1) hydrogen bonded, (2) vdW 
bonded, and (3) nonbonded in 2D.
The non-interacting state. In this third possible 
state, the test water has no interaction with its clockwise 
neighbor
 (19)
The same way as for the other two states, the isother-
mal-isobaric partition function is obtained by integrating 
over translational degrees of freedom
 (20)
510 Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
where v0 is the volume available to the test molecule in this 
state and is calculated as the van der Waals gas approxi-
mation
 (21)
Where vd = vs for MB water. Integrating of the parti-
tion function gives 
 (22)
These three expressions, Equations (14), (18) and 
(22), give the isobaric-isothermal ensemble Boltzmann 
weights of the three possible states of each water mole-
cule. We assume a mean-field attractive energy, –Na/v,52,53 
among cages, where a is the van der Waals dispersion pa-
rameter (0.02, here) and v is the average molar volume, 
which we will define below. The partition function for a 
single full hexagon of 6 waters would be given by
 (23)
Here, we treat hexagons a little differently instead. 
We define cooperative HB state or solid state that involves 
a higher degree of HB cooperativity than the hydrogen 
bonding that is just formed pairwise among nearest neigh-
bor waters in the liquid state. So, the partition function for 
each hexagon will be given by
 (24)
where δ = e–βec is the Boltzmann factor for the coopera-
tivity energy, ec, that applies only when 6 water molecules 
all connected into a full hexagonal cage. The terms on the 
right-side of this expression simply replace the statistical 
weight for each weakly hydrogen-bonded full hexagonal 
cage with the statistical weight for a cooperative strong-
ly hydrogen-bonded hexagonal cage. Δs is the Boltzmann 
factor for a cooperative HB or solid state. It differs from 
ΔHB only in that the former uses the hexagonal cage vol-
ume per molecule, vs, while the HB state uses the liquid 
water hydrogen bonding volume per molecule, vHB. Now 
we combine the Boltzmann factors for the individual wa-
ter molecules to get the partition function Q for the whole 
system of N particles,
 (25)
The factor N/6 accounts for the 3 possible interaction 
sites per water molecule and corrects for double counting 
the hydrogen bonds. We compute the populations of the 
states i = 1 (HB), 2 (LJ), 3(o), 4(solid) using
 (26)
The chemical potential is given by
 (27)
The molar volume is
 (28)
and all the other thermodynamic properties below are ob-
tained as described previously.52,53 For all the model cal-
culations, we used the parameters from potential function 
eHB, eLJ, rHB and sLJ. The parameter k = 10 was determined 
from the shape of the MB potential while ec = 0.03 was 
determined empirically.
3. 2. 3D MB Model
Here, we will point out only the differences between 
the theory in 3D in comparison to 2D.38 We consider a 
system of 3D MB model water molecules, modeled as 
three-dimensional spheres, and suppose that the structure 
of the liquid state of 3D model water is a perturbation from 
an underlying hexagonal (ice) lattice; (See Figure 5). Each 
molecule can be in one of the three possible orientational 
states like in 2D. These states are graphically presented in 
Fig. 6. 
Figure 5: Lattice of the model showing both the hexagon of the ice-
like structure and a pair interaction used for bookkeeping to avoid 
triple counting. Presented is only one layer.
The hydrogen-bonded state. If the test water mole-
cule points one of its four hydrogen bonding arms at an 
angle θ to within π/3 of the center of its clockwise neigh-
bor water, it forms a hydrogen bond. This is equivalent to 
about one fourth of the full solid angle. The energy of in-
teraction of the test water is
 (29)
k is the angular spring constant that describes the weak-
ening of the hydrogen bond as it becomes increasingly off 
angle. To compute the isothermal-isobaric partition func-
511Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
tion, ΔHB, of this HB state, we integrate this Boltzmann 
factor over all the allowable angles and over all the allowa-
ble separations of the test molecule relative to its clockwise 
neighbor,
 (30)
Where c(T) is here the 3D version of the kinetic 
energy contribution to the partition function. òòò dxfydz 
represents the volume over which the second molecule has 
translational freedom to form a hydrogen bond with the 
first water and is equal to effective volume veffHB = 0.242. 
The double integral òò dafy sums the orientations over 
which the test molecule has orientational freedom and 
is equal to 4π2. The volume vHB of the hydrogen-bonded 
state we determine similarly as for the 2D model. For the 
perfect hexagon crystal representing low-pressure ice, the 
volume of the solid is
 (31)
We estimate volume vHB as perturbation of this state as
 (32)
where xv = 1.12 is chosen empirically because density of 
the liquid state at room temperature is about 12% more 
dense than ice. Using these definitions and performing the 
integration in Equation (30) gives
 (33)
Figure 6: The three model states: (1) hydrogen bonded, (2) vdW 
bonded, and (3) nonbonded in 3D.
The van der Waals (vdW) state. Here, the test wa-
ter molecule forms only a van der Waals contact with its 
clockwise neighboring water. The water molecule has an 
energy
 (34)
The isothermal-isobaric partition function, ΔLJ is
 (35)
The triple integral òòò dxfydz represents the transla-
tion volume over which the second molecule forms a van 
der Waals contact with the first water and is equal to ef-
fective volume veffLJ = 0.086. The integrals over angles are 
equal to 8π2. The volume vLJ of this state is approximated 
as a volume of cubic close-packed crystal where the closest 
molecules are at distance σLJ√
6     -2
 (36)
we also tried other symmetries, but the results did not 
change much. Integrating of the partition function gives 
 (37)
The non-interacting state. In this third possible 
state, the test water has no interaction with its clockwise 
neighbor and the isothermal-isobaric partition function is 
equal to
 (38)
Here, we also assume a mean-field attractive energy, 
–Na/v,52,53 among cages, where a is the van der Waals dis-
persion parameter (0.045, here). The partition function for 
a single full hexagon of 6 waters and other properties are 
calculated in the same way as in 2D. For all the model cal-
culations, we used the parameters from potential function 
eHB, eLJ, rHB and sLJ. The parameter k = 80 was determined 
from the shape of the 3D MB potential while ec = 0.18 was 
determined empirically.
3.3. The Real Water – CageWater
Here we made slight modification in comparison with 
3D MB.54 Two water molecules can interact through a hydro-
gen bond (which depends on their relative orientations), inter-
act through a contact (which is orientation independent and 
occurs when they are close in space and no HB is present), or 
be noninteracting (when they are far apart, as in van der Waals 
gas). Hydrogen bonds are further parsed into two types: an 
HB can occur between 2 adjacent waters that have no higher 
order structure or can occur within a 12-water hexagonal unit 
cell (cage) forming 15 HBs. Parameters needed for calcula-
tions were obtained by getting good agreement with tempera-
512 Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
ture dependence of density at normal pressure and of boiling 
point position and are presented in Table 1.
4. Results
We present our results below in dimensionless or 
reduced units for both MB models, normalized to the 
strength of the optimal hydrogen bond eHB and hydrogen 
bond separation rHB  
for 2D MB model and  for 3D MB 
model).
4. 1. 2D MB Model
Analytical theory has additional approximations 
compared to computer simulations, which is why we first 
Table 1: To obtain the parameters, the intrinsic HB energy and HB distance are fixed while all 
other parameters were optimized.
Parameter Value Description
eHB 4.106 kcal/mol intrinsic HB energy between two molecules
rHB 2.767 Å intrinsic HB distance between two molecules
k 225.83 kcal/mol angle flexibility of HB
eLJ 0.8212 kcal/mol intrinsic LJ energy between two molecules
sLJ 3.293 Å intrinsic LJ distance between two molecules
vd 16.85 Å3 hard core of water molecule
ec –0.252 kcal/mol correlation energy per bond in 12-mer
xv 1.133 ratio between volumes of strong and weak HB states
veffHB 42369.9 Å3 effective volume of HB state
veffs 48089.8 Å3 effective volume of s state
veffLJ 74147.3 Å3  effective volume of LJ state
Figure 7: Temperature dependence of the density (ρ*), heat capacity (cp*), thermal expansion coefficient (α*), and isothermal compressibility (κ*) for 
2D MB water for different pressures. Results from the theory are plotted with lines and from the computer simulations40 by points.
513Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
checked the quality of the predictions of the analytical 
theory. We calculated the temperature dependence of the 
density (ρ*), heat capacity (cp*), thermal expansion coef-
ficient (α*), and isothermal compressibility (κ*) for dif-
ferent pressures. For a 2D MB model it was previously 
shown that the Mercedes-Benz water qualitatively cor-
rectly reproduces the anomalies of water27,31,32 for these 
quantities. In Fig. 7 a comparison of predictions of the 
present theory (lines) for the density, the thermal expan-
sion coefficient, the isothermal compressibility, and the 
heat capacity vs temperature to NPT Monte Carlo sim-
ulations (symbols) of the 2D MB model with the same 
parameters is shown. The calculations of the theory were 
performed at reduced pressure of 0.08, 0.12, 0.19 and 
0.32. The theory is in good general agreement with the 
simulations, including the density maximum (minima 
in molar volume). The thermal expansion coefficient is 
negative at low temperatures, which is consistent with 
computer simulations and with experiments for water. 
The Monte Carlo simulations of MB water do not show 
an experimentally observed minimum in the isothermal 
compressibility versus temperature. On the other hand, 
the present theory predicts a minimum for higher pres-
sures. At low temperatures, our present model shows a 
drop in heat capacity as the temperature is reduced. Be-
ing satisfied with the prediction of the model, we calcu-
lated non crystalline part of the phase diagram, shown in 
Fig. 8. The 2D MB model exhibits two critical points: the 
liquid-gas critical point (C1) at temperature T* = 0.118 
and pressure p* = 0.00035 which is slightly lower than ob-
tained from computer simulations,40 and the liquid-liq-
uid critical point (C2) at temperature 0.0212 and pressure 
0.42. There exists also a region of pressures between both 
critical points where we have only one fluid phase, at 
higher pressures we have two liquid phases, and at lower 
pressures the liquid and the gas phases. 
Figure 9: Temperature dependence of the diffusion (D*), viscosity (h*), thermal conductivity (l*) and thermal diffusivity (ld*) for 2D MB water for 
different pressures calculated by the theory. 
Figure 8: Phase diagram of the noncrystalline phases of water. Red 
line is liquid-liquid and blue line liquid-gas coexistence line.
514 Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
We have also developed the theory for dynamical 
properties. The diffusion processes occur in fluid or gas 
whenever a property is transported in a manner resem-
bling a random walk. If we assume that the water mole-
cules are doing random walk, we can say, for our 2D mole-
cules in each state, that the diffusion is proportional to the 
step size, l, and a step frequency, n. The step frequency is 
proportional to the Boltzmann factor for the change of en-
ergy from bonded to free state. This means that the energy 
of interaction is negative of bonding energy of molecule. 
We assumed that the step size is equal to the characteristic 
length of interaction in each state (lHB = ls = rHB for HB 
and s state, lLJ = sLJ for LJ state and for 0 state for 0 state the 
average distance between molecules in this state l0 = √
–
v0). 
For our model, we calculated the diffusion constant as a 
sum of all states of individual contributions 
 (39)
Where Di = λI ni for HB, s, LJ and 0 state of water. The 
step frequency is equal to Boltzmann factor of negative av-
erage bonding energy (ui) 
 (40)
To model viscosity, h, we start from Stokes-Einstein 
relation between viscosity and diffusion coefficient, D,
 (41)
We can express viscosity from this equation as
 (42)
We see that we can calculate viscosity from diffusion 
coefficient, temperature and diameter, d, of particle. In our 
case, we use averaged particle diameter which we calculat-
ed as a sum over all possible states of water
 (43)
For water molecules in states HB, LJ and 0 we used 
diameter of molecule equal to rHB while for state s waters 
form hexagons and we use the diameter of water in hexa-
gon state as equal to 2rHB. Next, we calculated the speed of 
sound cs as
 (44)
and thermal conductivity l using modified Eyring’s for-
mula as
 (45)
and thermal diffusivity ld as
Figure 10: Pressure dependence of the diffusion (D*), viscosity (h*), thermal conductivity (l*) and thermal diffusivity (ld*) for 2D MB water for 
different temperatures calculated by the theory. 
515Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
 (46)
In Fig. 9 and 10, we have plotted temperature and 
pressure dependence of the dynamical properties (diffu-
sivity, viscosity, thermal conductivity, and thermal diffu-
sivity). All the quantities have similar anomalous non-
monotonic behavior as for experimental water.
4. 2. 3D MB Model
As for 2D MB we first checked the quality of the predic-
tions of the analytical theory also for 3D MB. We also calcu-
lated the temperature dependence of the molar volume, heat 
capacity, isothermal compressibility, and thermal expansion 
coefficient for different pressures. For a 3D MB model, it 
was previously shown that the Mercedes-Benz water quali-
tatively correctly reproduces the anomalies of water49,50,51 for 
these quantities. In Fig. 11, a comparison of predictions of 
the present theory (lines) for the molar volume (V*/N), heat 
capacity (cp*), thermal expansion coefficient (α*), and iso-
thermal compressibility (κ*), vs temperature to NPT Monte 
Carlo simulations (symbols) of the 3D MB model with the 
same parameters is shown. The calculations of the theory 
were performed at reduced pressure of 0.12 and 0.19. The 
theory is in good general agreement with the simulations, 
including the minima in molar volume (density maximum). 
The thermal expansion coefficient is negative at low temper-
atures, which is consistent with computer simulations and 
with experiments for water. The Monte Carlo simulations 
of MB water do not show an experimentally observed min-
imum in the isothermal compressibility versus temperature. 
On the other hand, the present theory predicts a minimum 
for higher pressures. At low temperatures, our present model 
shows a drop in heat capacity as the temperature is reduced. 
Being satisfied with the prediction of the model, we contin-
ued our research by calculating the density of 3D MB water as 
a function of temperature along isobars (up to p* = 0.25) and 
determine critical points of the model. Results are shown in 
Fig. 12. In this pressure range, upon increase of temperature 
density increases, reaches a maximum, and then decreases. 
We determined non crystalline part of the phase diagram, 
shown in Fig. 13. The 3D MB model exhibits two critical 
points: the liquid-gas critical point (C1) at temperature T*= 
0.117 and pressure p* = 0.0115, and the liquid-liquid critical 
point (C2) at temperature 0.0779 and pressure 0.167. There 
exists also a region of pressures between both critical points 
where we have only one fluid phase, at higher pressures we 
have two liquid phases, and at lower pressures the liquid and 
the gas phases. 
4. 3. The Real Water – CageWater
Here, we compare the measured properties over wa-
ter’s liquid range to those predicted by parametrization for 
Figure 11: Temperature dependence of the molar volume (V*/N), heat capacity (cp*), thermal expansion coefficient (α*), and isothermal compress-
ibility (κ*), for 3D MB water for pressures p* = 0.19 (orange) and p* = 0.12 (blue). Results from the theory are plotted with lines and from the com-
puter simulations by symbols.50,51
516 Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
experimental data, called CageWater.54 In Fig. 14a we have 
plotted experimental data8 and data by best practices wa-
ter simulation models TIP4P/2005,21 TIP3P,20 SPC,20 and 
mW56 of the four main thermal and volumetric properties 
of water: the density, the thermal expansion coefficient, 
the isothermal compressibility, and the heat capacity. The 
comparison to experiments shows that the present model 
gives equal or better agreement than the simulation mod-
els over the normal and supercooled liquid temperature 
range and does not have the fluctuation errors that simula-
tions have, but the theoretical model has more parameters 
The model allows us to parse the experimental observa-
bles into hydrogen-bonding, caging, van der Waals, and 
non-interacting molecular components. Water is known 
to have a high heat capacity (ability to absorb thermal 
energy upon heating) among liquids of similar molecular 
size. The main conclusions from Figure 14b. are the fol-
lowing. In the normal liquid range, the high heat capacity 
comes from the breaking of two types of bonds: pairwise 
H bonds and Lennard−Jones-like contacts. Heating hot 
water near the boiling point leads to lower density, as it 
would for any LJ fluid, because heating hot water chang-
es the contact interactions more than the H bonds. Figure 
14c shows the same bulk properties as in Figure 14a except 
now computed as a function of pressure, not temperature. 
As increasing pressure squeezes water to become more 
compact (density increases and compressibility decreases), 
it crumples the hexagonal water cages breaking them into 
component pieces that just have pairwise water−water hy-
drogen bonding with little change to LJ and noninteract-
ing water populations. Pressure decreases the heat capacity 
(bond-breaking capability) because although it melts out 
some cages it is also “freezing in” some pairwise H bonds. 
The thermal expansion coefficient increases with pressure 
because pressure melts out the rigid cages into fragment-
ed H-bond pairs, which can be more readily squeezed 
together by pressure. Figure 15 shows that CageWater ac-
curately reproduces the anomalous hallmark thermal and 
volumetric signatures of the LLPT, namely, the divergent 
increasing heat capacity and compressibility with lowered 
temperature. Moreover, this model gives the microscopic 
components of those observables. We find that the large 
diverging heat capacity is due to the water cages, which 
have dominant populations in cold and supercooled water. 
The heat capacity is the sum of two contributions for each 
state: the population of that state multiplied by the indi-
vidual heat capacity. We also find that the negative thermal 
expansion of supercooled water is dominated by the cage 
term. Heating supercooled water shrinks the average vol-
ume by melting the cages, which are voluminous, and con-
verts them to smaller H-bonded fragments, like breaking 
a glass jar into shards that pack more compactly. This same 
physics is reflected in the peak of the compressibility at the 
supercooling peak temperature. Our model indicates that 
the two liquids that are in equilibrium around −50 °C are 
cage structures and broken H-bonded pieces.
5. Conclusion and Future  
Perspectives
We developed an analytical theory of water and ap-
plied it to 2D MB, 3D MB water models and parametrized 
it for experimental data. We used it for explaining how the 
pVT properties of liquid water arise from water’s hydrogen 
bonding and contacts. The theory predicts volumetric and 
energetic properties rather well, for experimental data it is 
more accurate than explicit simulation models yet is much 
faster to compute. Its simplicity and predictive power come 
from representing water using only four factors in the par-
tition function, hydrogen bonds, Lennard-Jones contacts, 
noninteracting terms and cooperative cages, rather than as 
a more extensive density expansion. The analytical theory 
advances our understanding of water’s structure−proper-
ty relations in showing that water’s long-known 2-density 
behavior is encoded in relatively infrequent cages which 
melt out strongly with temperature and pressure. This un-
Figure 12: Temperature dependence of the density for various pres-
sures (blue solid lines), high-density liquid–low-density liquid co-
existence line (red dashed line), liquid-gas density coexistence line 
(green dashed line), and maximum densities (pink dashed line).
Figure 13: Phase diagram of the noncrystalline phases of water. 
Blue solid line is liquid-liquid and orange dashed line liquid-gas co-
existence line.
517Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
Figure 14: Experimental liquid water properties (yellow full triangles) vs model predictions (dark blue solid line) and computer simulations. (a) 
Temperature dependence of liquid water’s density, thermal expansion coefficient, isothermal compressibility, and heat capacity at 1 bar pressure. 
Experiments are from Eisenberg and Kauzmann.8 Computer simulations are from Abascal and Vega for TIP4P/200521 and Jorgensen and Jenson for 
TIP3P20 and SPC20, and mW model predictions from Molinero and Moore56. (b) Molecular constituents of water at different temperatures: HB 
(hydrogen-bonded waters), s (12-mer hexagons), and LJ (waters in contact but not hydrogen bonded). (c) Pressure dependences of the same prop-
erties and their constituents at a temperature of 273 K.
518 Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
derstanding of water structure−property relations may aid 
in engineering filtration, osmosis, and desalination mate-
rials, in better solvation models for drugs and biomolecule 
actions, and for interpreting planetary geochemistry and 
hydrological cycles. The challenge of the current model is 
in calculating dielectric permittivity57,58 and solvation of 
polar molecules and ions. To be able to do this the model 
will have to be upgraded to version to include also charges 
on the water test particles, but this will add new param-
eters. Another challenge in theoretical modelling lies in 
developing the same kind of theory for other compounds 
like ionic liquids55, alcohols etc. We are expecting that all 
this can be done.
6. Acknowledgments
The financial support of the Slovenian Research 
Agency through Grant P1-0201 as well as to projects J7-
1816, J1-1708, N1-0186 and N2-0067 is acknowledged 
as well as National Institutes for Health RM1 award 
RM1GM135136.
Figure 15: Molecular components of supercooled water vary with temperature and pressure. Colored lines show the components HB (pairwise H 
bonded), s (12-mer cages), LJ (waters in pairwise contact), and 0 (waters separated and noncontacting). Dark blue line is the sum of all components. 
Pressure dependence is calculated at −35 °C and temperature dependence at 0.1 MPa (1 atm). Most definitive features are the strong variations of the 
balance of molecular components with T and p and how strongly the caging behaviors are opposed by the pairwise hydrogen-bonded waters.
519Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
7. References
 1  M. Chaplin, Nat. Rev. Mol. Cell Biol. 2006, 7, 861.
 DOI:10.1038/nrm2021
 2.  M. J. Tait, F. Franks, Nature 1971, 230, 91.
 DOI:10.1038/230091a0
 3.  P. Gallo, K. Amann-Winkel, C. A. Angell, M. A. Anisimov, F. 
Caupin, C. Chakravarty, E. Lascaris, T. Loerting, A. Z. Pana-
giotopoulos, J. Russo, J. A. Sellberg, H. E. Stanley, H. Tanaka, 
C. Vega, L. Xu and L. G. M. Pettersson, Chem. Rev. 2016, 116, 
7463.   DOI:10.1021/acs.chemrev.5b00750
 4.  E. Brini, C. J. Fennell, M. Fernandez-Serra, B. Hribar-Lee, M. 
Luksic, and K. A. Dill, Chem. Rev. 2017, 117, 12385.
 DOI:10.1021/acs.chemrev.7b00259
 5.  P. Žnidaršič-Plazl, Acta Chim. Slov. 2021, 68, 1.
 DOI:10.17344/acsi.2020.6488
 6.  J. Vladić, N. Nastić, T. Stanojković, Ž. Žižak, J. Čakarević, L. 
Popović and S. Vidovic, Acta Chim. Slov. 2019, 66, 473.
 DOI:10.17344/acsi.2019.5011
 7.  F. Franks., Ed. Water, a Comprehensive Treatise, (Plenum 
Press, New York, 1972–1980) Vol. 1–7.
 8.  D. Eisenberg and W. Kauzmann, The structure and properties 
of water (Oxford University Press, Oxford, 1969).
 9.  W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, 
and M. L. Klein, J. Chem. Phys. 1983, 79, 926.
 DOI:10.1063/1.445869
10 . I. Nezbeda, J. Mol. Liq. 1997, 73/74, 317.
 DOI:10.1016/S0167-7322(97)00076-7
11 . B. Guillot, J. Mol. Liq. 2002, 101, 219.
 DOI:10.1016/S0167-7322(02)00094-6
12.  C. Vega, J. L. F. Abascal, M. M. Conde, and J. L. Aragones, 
Faraday Discuss. 2009, 141, 251.   DOI:10.1039/B805531A
13.  F. H. Stillinger, Science 1980, 209, 451.
 DOI:10.1126/science.209.4455.451
14.  C. Tanford, The hydrophobic effect: formation of micelles 
and biological membranes, 2nd ed. (Wiley, New York, 1980,.
15.  W. Blokzijl and J. B. F. N. Engberts, Angew. Chem. Int. Ed. 
Engl. 1993, 32, 1545.   DOI:10.1002/anie.199315451
16.  G. Robinson, S.–B. Zhu, S.Singh and M. Evans, Water in Bi-
ology, Chemistry and Physics: Experimental Overviews and 
Computational Methodologies (World Scientific, Singapore, 
1996).   DOI:10.1142/2923
17.  R. Schmidt, Monatshefte fu¨r Chemie 1993, 132, 1295.
18.  A. Ben-Naim, Biophys. Chem. 2003, 105, 183.
 DOI:10.1016/S0301-4622(03)00088-7
19.  L. R. Pratt, Annu. Rev. Phys. Chem. 2002, 53, 409.
20.  L. W. Jorgensen and C. Jenson, J. Comput. Chem. 1998, 19, 
1179.
21 . J. L. F. Abascal and C. Vega, J. Chem. Phys. 2005, 123, 234505.
 DOI:10.1063/1.2121687
22.  E. Lambros and F. Paesani, J. Chem. Phys. 2020, 153, 060901.
 DOI:10.1063/5.0017590
23.  K. B. Lipkowitz, D. B. Boyd, S. J. Smith, and B. T. Sutcliffe, 
Rev. Comput. Chem. 2007, 10, 271.
24.  E. J. Baerends and O. V. Gritsenko, J. Phys. Chem. 1997, A 
101, 5383.   DOI:10.1021/jp9703768
24.  D. van der Spoel, P. J. van Maaren, and H. J. C. Berendsen, J. 
Chem. Phys. 1998, 108, 10220.   DOI:10.1063/1.476482
25.  H. W. Horn, W. C. Swope, J. W. Pitera, J. D. Madura, T. J. Dick, 
G. L. Hura, and T. Head- Gordon, J. Chem. Phys. 2004, 120, 
9665.   DOI:10.1063/1.1683075
26.  J. L.F. Abascal and C. Vega, J. Chem. Phys. 2005, 123, 234505.
 DOI:10.1063/1.2121687
27.  K. A. T. Silverstein, A. D. J. Haymet and K. A. Dill, J. Am. 
Chem. Soc. 1998, 120, 3166.   DOI:10.1021/ja973029k
28 . A. Ben–Naim, J. Chem. Phys. 1971, 54, 3682.
 DOI:10.1063/1.1675414
29 . A. Ben–Naim, Mol. Phys. 1972, 24, 705.
 DOI:10.1080/00268977200101851
30 . C. L. Dias, Phys. Rev. Lett. 2012, 109, 048104.
 DOI:10.1103/PhysRevLett.109.048104
31.  T. Urbic, V. Vlachy, Yu. V. Kalyuzhnyi, N. T. Southall and K. A. 
Dill, J. Chem. Phys. 2000, 112, 2843.   DOI:10.1063/1.480928
32.  T. Urbic, V. Vlachy, Yu. V. Kalyuzhnyi, N. T. Southall and K. A. 
Dill, J. Chem. Phys. 2002, 116, 723.   DOI:10.1063/1.1427307
33.  T. Urbic, V. Vlachy, Yu. V. Kalyuzhnyi and K. A. Dill, J. Chem. 
Phys. 2003, 118, 5516.   DOI:10.1063/1.1556754
34.  T. Urbic, V. Vlachy, O. Pizio, K. A. Dill, J. Mol. Liq. 2004, 112, 
71.    DOI:10.1016/j.molliq.2003.12.001 
35.  T. Urbic, V. Vlachy, Yu. V. Kalyuzhnyi and K. A. Dill, J. Chem. 
Phys. 2007, 127, 174511.   DOI:10.1063/1.2784124
36.  T. Urbic and M. F. Holovko, J. Chem. Phys. 2011, 135, 134706.
 DOI:10.1063/1.3644934
37.  T. Urbic and K. A. Dill, J. Chem. Phys. 2010, 132, 224507.
 DOI:10.1063/1.3454193
38 . T. Urbic, Phys. Rev. E 2012, 85, 061503.
 DOI:10.1103/PhysRevE.85.061503
39 . T. Urbic, Phys. Rev. E 2016, 94, 042126.
 DOI:10.1103/PhysRevE.94.042126
40 . T. Urbic, Phys. Rev. E 2017, 96, 032122.
 DOI:10.1103/PhysRevE.96.032101
41.  T. Urbic, V. Vlachy, O. Pizio, K. A. Dill, J. Mol. Liq. 2004, 112, 
7180.   DOI:10.1016/j.molliq.2003.12.001
42.  M. Kurtjak, T. Urbic, Mol. Phys. 112, 2014, 1132–1148.
 DOI:10.1080/00268976.2013.836608
43.  M. Kurtjak, T. Urbic, Mol. Phys. 2015, 113, 727–738.
 DOI:10.1080/00268976.2014.973919
44.  T. Urbic, V. Vlachy, K. A. Dill, J. Phys. Chem. 2006, B 110, 
49634970.   DOI:10.1021/jp055543f
45.  T. Urbic, M. F. Holovko, J. Chem. Phys. 2011, 135, 19.
 DOI:10.1063/1.4967807
46 . T. Urbic, T. Mohoric, J. Chem. Phys. 2017, 146 094505.
 DOI:10.1063/1.4977214
47 . P. Ogrin, T. Urbic, J. Mol. Liq. 2020, 114880.
 DOI:10.1016/j.molliq.2020.114880
48 . P. Ogrin, T. Urbic, J. Mol. Liq. 2021.
49.  C. L. Dias, T. Hynninen, T. Ala–Nissila, A. S. Foster and M. 
Karttunen, J. Chem. Phys. 2011, 134, 065106.
 DOI:10.1063/1.3537734 
50.  A. Bizjak, T. Urbic, V. Vlachy, and K. A. Dill, Acta Chim. Slov. 
2007, 54, 532.
51.  A. Bizjak, T. Urbic, V. Vlachy and K. A. Dill, J. Chem. Phys. 
520 Acta Chim. Slov. 2021, 68, 505–520
Urbic:  Analytical Modeling of Thermodynamic and Transport    ...
2009, 131, 194504.   DOI:10.1063/1.3259970
52.  T. M. Truskett and K. A. Dill, J. Chem. Phys. 2002, 117, 5101.
 DOI:10.1063/1.1505438
53.  T. M. Truskett and K. A. Dill, J. Phys. Chem. 2002, B 106, 
11829.   DOI:10.1021/jp021418h
54.  T. Urbic and K. A. Dill, J. Am. Chem. Soc. 2018, 140, 17106.
 DOI:10.1021/jacs.8b08856
55.  M. Bešter-Rogač, Acta Chim. Slov. 2009, 67, 1.
 DOI:10.17344/acsi.2020.5870
56.  V. Molinero and E. B. Moore, J. Phys. Chem. 2009, B 113, 
4008.   DOI:10.1021/jp805227c
57.  A. V. Dubtsov, S. V. Pasechnik, D. V. Shmeliova, A. Sh. Said-
gaziev, E. Gongadze, A. Iglič, S. Kralj, Electrostatic Effects, Soft 
Matter, 2018, 14, 9619–9630.   DOI:10.1039/C8SM01529E
58.  A. Iglič, E. Gongadze, V. Kralj-Iglič, Acta Chim. Slov. 2019, 
66, 534.   DOI:10.17344/acsi.2019.5495
Except when otherwise noted, articles in this journal are published under the terms and conditions of the  
Creative Commons Attribution 4.0 International License
Povzetek
Struktura in lastnosti biomolekul, kot so proteini, nukleinske kisline in membrane, so odvisne od vode. Voda je zelo 
pomembna tudi v industriji. Na splošno je nenavadna snov z več kot 70 anomalnimi lastnostmi. Zaradi teoretičnega 
in računalniškega modeliranja vode jo vse bolje razumemo. Obstajajo različni tipi modelov vode. Prvi so kompleksi, 
ki upoštevajo veliko podrobnosti, in njihova računska zahtevnost narašča z detelji. Drugi pa so enostavni, ki se osre-
dotočajo na razlago osnovnih zakonitosti kot so termodinamika, fazni diagram. Ti modeli so računsko manj potratni. 
Enostavni modeli omogočajo boljše razumevanje vode na način, ki dopolnjuje kompleksne modele. Tu predstavljamo 
analitično modeliranje lastnosti vode na različnih nivojih, dve- in tridimenzionalnega Mercedes Benz modela vode ter 
eksperimentalne vode.