IOP PUBLISHING MODELLING AND SIMULATION IN MATERIALS SCIENCE AND ENGINEERING

Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 (16pp) doi:10.1088/0965-0393/18/2/025007

A parallel reaction-transport model applied to cement

hydration and microstructure development

Jeffrey W Bullard1, Edith Enjolras2, William L George2,

Steven G Satterfield2 and Judith E Terrill2

1 Materials and Construction Research Division, National Institute of Standards and Technology,

Gaithersburg, MD, USA

2 Mathematical and Computational Sciences Division, National Institute of Standards and

Technology, Gaithersburg, MD, USA

E-mail: jeffrey.bullard@nist.gov

Received 6 July 2009, in final form 29 October 2009

Published 18 January 2010

Online at stacks.iop.org/MSMSE/18/025007

Abstract

A recently described stochastic reaction-transport model on three-dimensional

lattices is parallelized and is used to simulate the time-dependent structural

and chemical evolution in multicomponent reactive systems. The model,

called HydratiCA, uses probabilistic rules to simulate the kinetics of diffusion,

homogeneous reactions and heterogeneous phenomena such as solid nucleation,

growth and dissolution in complex three-dimensional systems. The algorithms

require information only from each lattice site and its immediate neighbors,

and this localization enables the parallelized model to exhibit near-linear

scaling up to several hundred processors. Although applicable to a wide range

of material systems, including sedimentary rock beds, reacting colloids and

biochemical systems, validation is performed here on two minerals that are

commonly found in Portland cement paste, calcium hydroxide and ettringite, by

comparing their simulated dissolution or precipitation rates far from equilibrium

to standard rate equations, and also by comparing simulated equilibrium states

to thermodynamic calculations, as a function of temperature and pH. Finally,

we demonstrate how HydratiCA can be used to investigate microstructure

characteristics, such as spatial correlations between different condensed phases,

in more complex microstructures.

1. Introduction

Understanding and predicting chemical and microstructural evolution in cement paste is a longstanding objective in concrete technology. The coupled dissolution and precipitation reactions

that occur when cement powder is mixed with water determine a number of properties that are

practically important, such as the setting time and the rate of heat release of concrete. Increased

0965-0393/10/025007+16$30.00 © 2010 IOP Publishing Ltd Printed in the UK 1
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

understanding of the mechanisms and rates of hydration reactions that occur in the first hours

after mixing, as well as the development of microstructure, could lead to advances in the design

and deployment of chemical admixtures for controlling these ‘early age’ properties. Better

ability to predict the early age properties of concrete could allow more informed decisions in

the field about how to transport and place the material in service.

Accurate prediction of early age properties of concrete with arbitrary mixture proportions

and chemical composition requires a model that can simulate the rates of the coupled

dissolution/precipitation reactions in a manner consistent with the development of the solid

microstructure. A previous publication describes the development of a stochastic, latticebased reaction-transport model that is a promising candidate for modeling the reactions and

microstructural development leading to setting in cementitious materials [1, 2]. That model,

called HydratiCA3, has been demonstrated to accurately simulate the diffusion of neutral or

charged species in dilute and semi-dilute solutions [1]. It has also been shown to simulate the

kinetics of simple unimolecular or bimolecular reactions, both homogeneously in solution and

heterogeneously at solid–liquid interfaces [2]. Recently, HydratiCA has been used to clarify

the mechanisms of the hydration of tricalcium silicate [3], a major component of Portland

cement. In contrast to earlier microstructure models developed specifically for cement paste

hydration, such as CEMHYD3D [4] and μIC [5], the model formulation enables simulation

of realistic chemical kinetics and thermodynamics and provides an intimate coupling to the

development of microstructure.

Earlier studies using HydratiCA invoked several simplifying assumptions about the system

for the sake of computational expediency. For example, to limit the number of coupled reactions

that were simulated, those earlier studies intentionally neglected the formation of ion complexes

in solution. In addition, past simulations necessarily have been confined to very small systems,

no more than 1000 lattice sites with a lattice spacing of 1μm, which is enough to model a

single cement particle in solution. These small sizes are justifiable for investigating purely

chemical kinetic effects. However, in cementitious materials and other complex porous media,

the representative volume elements required to compute microstructure properties are much

larger, typically by a factor of 100 or 1000. In such systems, the computational demands

are so high for HydratiCA that it cannot be validated properly unless it is modified to run in

parallel on multiple processors. Recently, this modification has been made, and we now use it

to examine larger-scale systems than have been possible before.

In the next section, the model and its parallelization are described briefly. In section 3,

we apply the model to several phenomena that are important in cement chemistry, including

the influence of temperature and pH on the dissolution kinetics and equilibrium of portlandite,

and the dissolution and equilibration of ettringite—a very sparingly soluble mineral in water.

Finally, we report simulation results on the microstructural development of an idealized model

of Portland cement paste, and the development of spatial correlations between growing and

dissolving condensed phases.

2. Modeling

The algorithms used in HydratiCA for modeling diffusion and chemical reactions are based

on fundamental kinetic cellular automaton principles [6]. The details and basic verification

of the algorithms are given elsewhere [1, 2]. The material microstructure is discretized on

a regular cubic lattice having a lattice spacing of λ. The initial cement particle and water

microstructure is mapped onto this lattice by assigning a phase (e.g. Ca(OH)2 or water) to each

3 The name HydratiCA is a pseudo-acronym for Hydration by Cellular Automata.

2
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

lattice site. These materials are themselves finely discretized into quanta of concentration

called cells; the number of cells of a given material at a particular lattice point determines its

local concentration.

Chemical changes and structural development are simulated by iterating over small time

increments, each of which is decomposed into independent transport and reaction steps.

Diffusion is modeled by allowing each cell at a lattice site to execute a random walk to a

neighboring site. The probability pt of the walk depends on the effective diffusivity D of the

mobile species at the site and the length of the time increment τ being considered [1],

pt = τD

λ2 . (1)

Similarly, probabilistic rules are formulated to simulate chemical reactions at a lattice site. The

probability of the reaction occurring depends on the reaction rate constant and on the number

of cells Nα of each reactant α that participates in the reaction. For example, the probability of

occurrence of the ith homogeneous reaction is given by [6]

p(i)

rx = kξ (



α ν(i)

α )−1

τ



α

max



0,

ν(i)

α

m=1

Nα − m + 1



 , (2)

where ξ is a constant model parameter that relates Nα to the molar concentration of species α

and ν(i)

α is the molar stoichiometric coefficient of the reactant α participating in the reaction.

Equation (2) strictly applies only to homogeneous reactions. The same kind of equation applies

for heterogeneous reactions and nucleation events, although the length scaling is somewhat

different [2]. The reaction is allowed if p(i)

rx exceeds a random number q ∈ [0, 1] drawn from a

uniform distribution. If the reaction occurs at a lattice site, the number of cells of each reactant

(product) is decremented (incremented) by the number required by the molar stoichiometric

coefficients of the reaction. Reactions that can proceed at appreciable rates in both forward

and reverse directions are modeled as two separate one-way reactions. Both equations (1)

and (2) have been shown to converge, in the limit τ , λ → 0, to the standard rate equations for

diffusion and homogeneous reactions, respectively [6].

Besides the molar stoichiometric coefficients for the reactants and products, each reaction

is characterized by its absolute rate constant in either the forward or the reverse direction,

the solubility product or equilibrium constant at a reference temperature (taken to be 298 K

unless stated otherwise), the activation enthalpy for either the forward or reverse reaction, the

enthalpy of reaction and nucleation energy barriers (for heterogeneous reactions only). These

parameters are sufficient to capture the rate of reaction, the temperature dependence of the rate

and even the temperature dependence of the equilibrium state for reversible reactions [2]. In

addition to these reaction parameters, several properties of each of the constituent substances

are also required. For simple condensed phases such as water and stoichiometric solids, the

density, molar volume and diffusion formation factor must be supplied. For mobile ionic solute

species, one must provide the electrical charge, the diffusion coefficient at infinite dilution and

the Kielland ion-size parameters needed to calculate the activity coefficients using the extended

Debye–Huckel equation. Values for many of these properties can be found in textbooks or  ̈

other reference materials [7–10]. Unless stated otherwise, all equilibrium constants, Kielland

ion-size parameters and other thermodynamic data used in this paper are taken from the

Nagra/PSI chemical thermodynamic database [11]. Finally, for non-stoichiometric solids,

the model simulates compositional and structural variability by microscopic coprecipitation

of two stoichiometric end member phases that span the desired range of compositions and

physical properties, as described more fully in [3]. Each end member is assigned values for

all the properties just listed, and then the relative rates of formation of each end member at a

3
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

lattice site determines the local composition of the phase. An approach similar to this is also

used in the thermodynamic modeling of non-stoichiometric solids [12].

2.1. Parallelization

Each lattice site in HydratiCA is occupied by numbers of cells of differing types, each of which

represents a discrete concentration of a particular material component or phase. The algorithms

used to compute changes in a site’s cell occupation numbers require information only from the

site and its nearest and next-nearest neighbor sites. This localization of the algorithms leads

to a natural parallel decomposition of the simulation, by which the lattice is partitioned into

sublattices, each of which is assigned to a separate computational process. Each process uses

one CPU processor. A process owns each of the lattice sites of the sublattice it is assigned, and

the owner of a lattice site is responsible for all of the computations needed to update that lattice

site during each time step of the simulation. Therefore, lattice sites that are on the surface of

a sublattice must gather data from neighboring lattice sites that are owned by other processes.

This updating is accomplished by surrounding each sublattice with a layer of ‘ghost’ lattice

sites that are images of the actual lattice sites owned by adjacent processes. In the simulations

reported here, the lattice was decomposed into orthorhombic sublattices that were as close to

the same size and as close to cubic as possible for the lattice shape and number of processes,

in an effort to minimize the number of ghost lattice sites. Inter-process communication is

handled using Message Passing Interface (MPI), the dominant communications protocol used

in high-performance parallel computing [13].

If the communications needed to update the ghost lattice sites in each time step required

no computational time (i.e. ideal communications), then simulations would scale perfectly up

to the point at which each process owned a single lattice site. Because ideal communications

are unattainable, parallel execution will be efficient only if each process owns enough lattice

sites to make the overhead of communication time a small percentage of the overall simulation

time.

HydratiCA’s scaling behavior was examined using up to 448 processes on the

Columbia supercomputer at the National Aeronautics and Space Administration Advanced

Supercomputing Division (NASA/NAS). Figure 1(a) compares the multiplicative speed-up in

identical simulations as a function of the number of processes. The simulations involved

17 types of cells with 12 coupled nonlinear reactions on a computational domain with

100×100×100 lattice sites. Further details of the simulations are discussed in section 3.3. In

figure 1(a), all timings are normalized by the timing for the simulation run on one processor.

The figure indicates that HydratiCA can run effectively on this 100 × 100 × 100 system

using up to at least 448 processes. However, with increasing numbers of processes, the actual

speed-up diverges from the ideal linear scaling, and little or no performance is gained beyond

400 processes. This can be seen more clearly by plotting the same data in terms of the

process efficiency, as shown in figure 1(b). Efficiency is defined here as the percentage of the

simulation time that each process devotes to the actual simulation as opposed to inter-process

communication.

It should be noted that the performance indicated in figure 1 is dependent on the number of

lattice sites in the system and also on the chemical complexity, that is the number of independent

chemical components and chemical reactions that must be considered. If either the number of

lattice points or the chemical complexity increases, the performance will scale better to higher

numbers of processes than indicated in the figure, which was generated using a relatively

simple chemical system.

4
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

100

(a) (b)

200 300 400 500

Number of Processors, N

100

200

300

400

500

Speed-up

0 100 200 300 400 500

Number of Processors, N

30

40

50

60

70

80

90

100

Efficiency (%)

Figure 1. (a) Parallel speed-up, computed as the time required to complete the computation on

one processor divided by the time required to complete the computation on N processors. The

simulations were run on a 100 ×100 × 100 lattice of ettringite particles in aqueous solution. The

dashed line is the idealized linear speed-up behavior. (b) Process efficiency for the same simulations

as in (a). Efficiency with N processors is computed as the speed-up with N processes divided by

N. This is a measure of the overhead incurred as a result of managing the parallelization of the

algorithm, which includes, but is not limited to, inter-processor communication costs.

3. Results and discussion

3.1. Portlandite in water: effect of temperature

To validate the model implementation on a relatively simple aqueous mineral system, we

assume a one-step reaction for the reversible dissolution of portlandite, Ca(OH)2, in water,

along with the speciation of calcium hydroxide in solution according to

Ca(OH)2

k1

 Ca2+ + 2OH−, log Ksp = −5.20, (3)

Ca2+ + 2OH− k2

 CaOH+, log Keq = 1.22. (4)

The values of the solubility product Ksp and complexation constant Keq are those reported at

298 K in [11]. We assume a dissociation rate constant of k1 = 7.2 μmol(m2 s)−1, inferred

from published studies of the kinetics of portlandite crystal growth [14]. Also based on that

data, we assume an activation enthalpy for portlandite dissolution of H∗

1 = 57.5 kJ mol−1.

The enthalpy of dissociation is H1 = −17.06 kJ mol−1. A value for the enthalpy of the

complexation reaction could not be found in any of the thermodynamic databases we searched,

perhaps because the heat capacity and its temperature dependence for CaOH+ have not been

measured. Here we assume a value of H2 = −2.16 kJ mol−1, which is comparable to the

enthalpies of similar speciation reactions [12]. Similarly, we have not found a published value

of the rate constant k2 for formation of the CaOH+ complex; however, complexation reactions

in solution generally equilibrate rapidly compared with mineral dissociation and growth [15].

In the absence of any data, we chose to set the rate constant considerably greater than k1 to

promote near-equilibrium conditions in solution. However, if the rate constant is set too high,

smaller time steps are required to maintain stability of the algorithm embodied in equation (2).

Through a series of test simulations, we found that a rate constant, k2 = 16.6 mol(m3 s)−1, was

great enough to maintain near-equilibrium conditions in solution without requiring significantly

smaller time steps.

5
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

Figure 2. Suspension of portlandite particles in water. The cubic lattice is 100μm on each side,

with a lattice spacing of 1μm.

(This figure is in colour only in the electronic version)

The system, shown in figure 2, is a cubic volume 100μm on each side, with a lattice spacing

of 1μm, containing a suspension of portlandite particles in pure water. The initial solid volume

fraction is 0.237, and the initial specific surface area of the solid is 1.43μm2 μm−3, which

corresponds to a mean equivalent spherical diameter of 4.2μm.

Figure 3 shows the predicted total calcium concentration as a function of time at 293, 298

and 303 K. The data follow the expected qualitative trends. That is, the initial rate of dissolution

increases with temperature, as expected for any thermally activated process. In addition,

the total calcium concentration near equilibrium decreases with increasing temperature, also

as expected for any net exothermic dissolution process. Quantitatively, the initial rates of

dissolution, calculated from the limiting value of the slopes at zero time in figure 3, are shown

in table 1 and compared with the ‘target’ value that is inferred from the assumed values of k1

and H∗

1 . The calculated rate at 293 K exceeds the true rate by only 3.5%, and the discrepancy

decreases with increasing temperature, although at the two higher temperatures the calculated

initial rates are slightly lower than the expected values.

To more clearly depict the accuracy with which the model predicts the equilibrium

state of the system, the inset of figure 3 compares the near-equilibrium total calcium

concentration, [Ca2+] + [CaOH+] (solid), with the recommended values (dashed) published

by the International Union of Pure and Applied Chemistry [16]. The greatest discrepancy,

occurring at 293 K, is 0.9 mmol L−1 or 4%. We also have compared the total calcium

concentration, as well as [Ca2+] and [CaOH+] individually, with the values predicted by

the GEMS-PSI thermodynamic modeling software package [17], which uses the Nagra/PSI

thermodynamic database [11]. Again, the maximum discrepancy is less than 1 mmol L−1

(about 4%) and occurs at 293 K.

6
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

0 25 50 75 100 125 150 175 200

time (s)

0

5

10

15

20

25

[Ca2+] + [CaOH+] (mmol/L)

293 K

298 K

303 K

150 160 170 180 190

19

20

21

22

23

Figure 3. Prediction of the time dependence of the total calcium concentration in solution for

dissolution of portlandite in pure water at a constant temperature of 293, 298 or 303 K. Dashed lines

in the inset represent experimentally measured values of equilibrium total calcium concentration

reported in [16].

Table 1. Comparison of calculated and target initial dissolution rates of portlandite in water at 293,

298 and 303 K.

Calculated rate Target rate Difference

T (K) (μmol m−2 s−1) (μmol m−2 s−1) (%)

293 5.00 4.83 3.5

298 7.00 7.19 −2.6

303 10.50 10.51 −0.1

It is important to note that different thermodynamic modeling software packages and

databases may not agree on the equilibrium state of a system, even for one as simple as

portlandite in water. For example, the thermodynamic calculations just described were

repeated with the PHREEQC thermodynamic modeling software package using the LLNL and

WATEQ4 databases, hereafter referred to collectively as P/LW. At 298 K, P/LW agrees to within

1 mmol L−1 with the cited experimental data, with the GEMS software package calculation

and with the HydratiCA calculation of total calcium concentration. However, at 293 K, P/LW

calculates a total calcium concentration of 29 mmol L−1, which is about 7 mmol L−1, or about

33%, higher than either the experimental value, the value calculated by GEMS, or the value

predicted by HydratiCA. Similarly, at 303 K, P/LW calculates a total calcium concentration of

14 mmol L−1, which is about 6 mmol L−1, or about 30%, lower than any of the other values.

Clearly, the temperature dependence predicted by P/LW near 298 K for this system, although in

the correct direction, is much too severe, and is most likely due to its approximation of activity

coefficients using the Davies equation instead of the extended Debye–Huckel equation as  ̈

implemented by GEMS and by HydratiCA.

3.2. Portlandite in water: effect of pH

In Portland cement pastes, the solution in the capillary pore space typically has a high pH (12

to 13), the value depending in part on the concentration of readily soluble alkali salts in the

7
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

0 25 50 75 100 125 150 175 200

time (s)

0

5

10

15

20

25

30

[Ca2+] + [CaOH+] (mmol/L)

[NaOH] = 0.8 mmol/L

[NaOH] = 8.7 mmol/L

[NaOH] = 43 mmol/L

Figure 4. Prediction of the time dependence of the total calcium concentration in solution for

dissolution of portlandite in solutions with different concentrations of NaOH. Dashed lines indicate

expected equilibrium total calcium concentrations calculated by GEMS-PSI thermodynamic

modeling software [17].

cement. Therefore, it is important for a model of cement hydration to adequately capture the

effect of solution pH on the dissolution kinetics and equilibrium solubility of minerals. As

a prototype example, we consider the behavior of portlandite in aqueous solutions of NaOH.

The same system is used as in the last section, except that the particles are now suspended

in an aqueous solution of NaOH at concentrations of either 0.87, 8.7 or 43 mmol L−1 at a

constant temperature of 298 K. In addition to reactions (3) and (4), we also need to consider

the formation of the NaOH0 complex in solution according to

Na+ + OH− k3

 NaOH0, log Keq = −0.18. (5)

Again, to ensure that this complexation reaction is sufficiently rapid to keep it from being

rate-controlling, we assume k3 = 16.6 mol(m3 s)−1, the same value used for the CaOH+

complexation reaction.

Figure 4 shows the time dependence of the total calcium concentration at three different

NaOH concentrations. The first noteworthy feature of the plot is the near-equilibrium

predictions of total calcium concentration, compared with expected values calculated by

GEMS-PSI (dashed lines). The simulation values exceed the expected values by about

0.001 mmol L−1 at all three concentrations. The reason for this small but systematic

discrepancy is not known, although it appears to be about the same magnitude and sign as

for portlandite in pure water from the last section.

As a check on the kinetic behavior, the instantaneous rate of dissolution at zero time is

expected to be exactly the same as in the absence of NaOH because [Ca2+], and therefore the

ionic activity product (IAP) for portlandite, is still zero initially. Indeed, the mean dissolution

rate over the first 0.5 s of dissolution is equal to the predicted value in the second column

of table 1 at 298 K. It is clear from figure 4, however, that the rates become unequal within

seconds, which also is expected because increased concentrations of OH− with increasing

concentrations of NaOH result in a higher IAP as soon as the concentrations of the other

species are nonzero.

8
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

3.3. Solubility of ettringite in water

Having validated the model implementation in terms of the influence of temperature and pH

variations on a simple mineral, we now test the model on a more chemically complex mineral.

For this purpose, we choose the mineral ettringite, for which the dissociation reaction may be

modeled as [12]

Ca6Al2(SO4)3(OH)12 · 26H2O

k4

 6Ca2+ + 2Al(OH)

−

4 + 3SO2−4

+ 4OH− + 26H2O, log Ksp = −45.09 (6)

which has a reaction enthalpy of H4 ≈ 204.5 kJ mol−1 at 298 K [12, 18]. Ettringite is an

important mineral in Portland cement both as a product of the reaction between tricalcium

aluminate and calcium sulfate, and also at later times where it participates in chemical

degradation processes [8]. It provides a good test of the model because its solubility product

is nearly 40 orders of magnitude lower than that of portlandite. The concentrations of calcium

and hydroxyl ions in solution are both an order of magnitude less than the equilibrium

concentrations for portlandite and, in addition, aluminate and sulfate ions are also present

in even lower concentrations. Thermodynamic modeling calculations indicate that more than

a dozen solute species should be present in solution at equilibrium with ettringite, six of

which—Ca2+, CaOH+, SO2−

4 , CaSO0

4, Al(OH)−

4 and OH−—should be present in greater than

trace concentrations. Therefore, we must account for a second speciation reaction, in addition

to reaction (4),

Ca2+ + SO2−

4

k5

 CaSO0

4, log Keq = 2.44. (7)

In the absence of reported rate data, we once again assume k5 = 16.6 mol(m3 s)−1 to prevent

this speciation reaction from being rate-controlling.

The system used is the same as shown in figure 2, except that the solid phase is now

identified as ettringite instead of portlandite. Table 2 shows the near-equilibrium concentrations

of solute components at 298 K and compares them with the concentrations predicted by the

GEMS-PSI thermodynamic modeling application [17] assuming the same solubility product.

It should be noted that solid Al(OH)3, gibbsite, is thermodynamically expected to precipitate

in small quantities. However, for the purposes of this calculation we have neglected the

formation of gibbsite both in HydratiCA and in the thermodynamic calculations. As the table

indicates, the predicted concentrations are well matched to the thermodynamic calculations,

with no more than 0.03 mmol L−1 absolute difference. The values given in table 2 also are

consistent with the range of solution compositions found in available experimental studies

[19–22]. Therefore, even for sparingly soluble minerals that may be in equilibrium with solute

components at concentrations as low as several μmol L−1, the model implementation provides

good calculations of equilibrium and kinetics.

3.4. Phase correlations in hydrating cement paste

By tracking the development of the microstructure of a 3D system, one could begin to

investigate spatial correlations between growing and dissolving phases. Some solute species

in cement, such as aluminate ions, are added to solution by dissolution of tricalcium

aluminate (C3A)4 and tetracalcium aluminoferrite (C4AF), but not by tricalcium silicate (C3S).

These aluminate ions are needed for the growth of minerals such as ettringite or calcium

monosulfoaluminate, but not for portlandite. Silicate ions are added to solution by dissolution

4 In this section, we adopt conventional cement chemistry notation by which single letters are used to denote simple

oxide components of cement minerals, e.g. C = CaO, A = Al2O3, F = Fe2O3, S = SiO2, H = H2O.

9
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

Table 2. Comparison of calculated and expected values of the solute concentrations in equilibrium

with ettringite at 298 K, neglecting the possibility of gibbsite formation.

Calculated Expected Difference

Component (mmol L−1) (mmol L−1) (mmol L−1)

Ca2+ 1.83 1.80 0.03

OH− 1.31 1.29 0.02

SO2−

4 0.85 0.84 0.01

Al(OH)

−

4 0.67 0.66 0.01

CaSO0

4 0.16 0.16 0.00

CaOH+ 0.04 0.03 0.01

Figure 5. Simulated suspension of tricalcium silicate, tricalcium aluminate, and calcium sulfate

dihydrate in water: (a) initial suspension microstructure, (b) microstructure after three hours

of hydration and (c) same microstructure with silicate phases subtracted to reveal the calcium

sulfate and calcium aluminate phases more clearly. The color scheme is light blue = C3S, light

brown = C–S–H, dark blue = portlandite, gray = C3A, yellow = gypsum, green = ettringite.

of C3S and C2S, but not from C3A or C4AF, and are needed for the growth of calcium silicate

hydrate gel (C–S–H) but not portlandite or ettringite. Therefore, a hydrating cement paste is

composed of a set of coupled reaction networks, similar to those being studied in biological

systems [23, 24], and the topology of these networks, along with the transport rates of the

solute species, may give rise to spatial correlations between different dissolving and growing

phases in cement paste.

As an example, we conclude with preliminary results for a system composed of C3S, C3A

and gypsum (CaSO4 · 2 H2O) suspended in water, shown in figure 5, that represents a highly

idealized Portland cement paste. The assumed reactions are given in table 3, along with the

required reaction properties, i.e. the forward rate constant, the equilibrium constant and the

enthalpy of reaction. The total solid volume fraction of the system is 0.32, which is somewhat

more dilute than typical Portland cement binders in concrete.

As the initial particles dissolve, the solution becomes supersaturated with respect to the

precipitation of solid hydration products, which form either on the surfaces of existing particles

or in the bulk solution. In particular, calcium silicate hydrate and portlandite are possible

hydration products of the dissolution of C3S (see reactions (1), (4)–(6) in table 3). Ettringite

and portlandite are possible hydration products formed by the dissolution of C3A and gypsum

(see reactions (2), (3), (6) and (7) in table 3). These two reaction networks are coupled through

Ca2+ and OH− ions in solution, which are common to both networks.

10
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

Table 3. Summary of assumed reactions and reaction parameters for hydration of idealized Portland

cement. k+ is the reaction rate constant in the forward direction (left to right) for each reaction.

k+ Hrx

Reaction (mol(m2 s)−1) log Keq (kJ/mol)

(1) C3S + 3H2O  3Ca2+ + H2SiO2−

4 + 4OH− 5.2 × 10−7a −17.0a −137.6b

(2) C3A + 6H2O  3Ca2+ + 2Al2O−

4 + 4OH− 2.0 × 10−6a −10.3b −248.3b

(3) CaSO4 · 2H2O  Ca2+ + SO2−

4 + 2H2O 2.0 × 10−5c −4.48b −0.46b

(4) C–S–H(I)  Ca2+ + H2SiO2−

4 + 3H2O 4.2 × 10−7a −7.52a 20b

(5) C–S–H(II)  2Ca2+ + H2SiO2−

4 + 2OH− + 3H2O 5.0 × 10−7a −12.96a 20b

(6) Ca(OH)2  Ca2+ + 2OH− 7.2 × 10−6d −5.2b −17.1b

(7) C3A · 3CaSO4 · 32H2O  6Ca2+3SO2−

4 + 2Al(OH)

−

4 + 4OH− + 26H2O 1.0 × 10−19a −45.1b 202.6b

(8) C3A · CaSO4 · 12H2O  4Ca2+SO2−

4 + 2Al(OH)

−

4 + 4OH− + 6H2O 1.0 × 10−14a −29.3b 45.6b

(9) C4AH19  4Ca2+ + 2Al(OH)

−

4 + 6OH− + 12H2O 1.0 × 10−9a −25.6b −0.53b

(10) Al(OH)3 + OH−  Al(OH)

−

4 2.22a −0.24b 22.3b

(11) CaOH+  Ca2+ + OH− 0.063a −1.2b −21.7b

(12) CaSO0

4  Ca2 + +SO2−

4 0.063a −2.1b −5.4b

(13) C3A + CaSO0

4  C3A + CaSO4(ads) 0.063a −2.1a 4.6a

a Assumed based on prior simulations. b Taken or estimated from thermodynamic tables in [9, 11]. c Inferred from published crystal growth data in [25]. d Inferred from published crystal growth data in [14].

Each reaction in table 3 results in the release or absorption of a specific quantity of heat

per mole of solid. A cursory examination of the enthalpies of reaction shown in the table

reveals that the major contributions to the heat evolved come from the dissolution of C3S and

C3A and the precipitation of ettringite (reaction 7). Therefore, the aggregate of the reactions

causes a well-defined heat signature that tracks the overall rate of reaction and that can be

measured by isothermal calorimetry. Isothermal calorimetry is commonly used in practice to

track the overall progress of the hydration reaction. The simulation tracks the incremental

progress of each reaction at each lattice site, so it automatically tracks the net enthalpy change

in the system and can be used to generate a simulated isothermal calorimetry curve, shown

in figure 6 as rate of heat released per gram of solid in the unhydrated system. The figure is

qualitatively typical of those observed for Portland cement pastes. Specifically, a large but brief

exothermic peak is observed within the first minutes of hydration due to rapid dissolution of

the cement particles after wetting. This initial reaction decelerates quickly and a period of slow

reaction, commonly called an induction period, can endure for minutes to hours depending

on the particle size distribution, temperature and specific composition of the cement. The

hydration reaction gradually accelerates as hydration products continue to grow. Eventually,

the hydration products impinge on each other and can partially or totally cover the dissolving

cement particles. The result is a second peak, much shorter and broader than the first peak, and

a gradual deceleration of the hydration process at later times. The maximum heat release rate

in this second peak, about 7 mW g−1, is consistent with values observed for Portland cement

pastes [26]. However, the length of the induction period shown in the figure is shorter than the

1 to 3 h that is commonly observed. Evidently, one or more of the reaction rate constants in

table 3 should be decreased, or the nucleation energy barriers need to be adjusted, for the length

of the simulated induction period to agree better with the typical behavior of Portland cement

pastes. More experimental research on simple systems is required to get better estimates of

some of the rate constants that have been assumed here.

Aside from the overall progress of the hydration reaction, we also are interested in the

development of microstructure in the cement, because an accurate microstructure enables

11
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

0 25 50 75 100 125 150

time (min)

0

5

10

15

20

Power (mW/g)

Figure 6. Calculated rate of heat release for the hydration of the C3S+C3A+ gypsum system

shown in figure 5.

subsequent analysis to estimate important macroscopic properties such as setting time, effective

elastic properties, and transport properties such as the permeability to fluids. Snapshots of the

microstructure at t = 0 h and at t = 3 h are shown in figure 5. In the rest of this section, we

will analyze the spatial distribution of the hydration products and their correlations with the

initial unhydrated phases.

A growing body of experimental [27–29] and modeling [3, 30] evidence strongly suggests

that C–S–H nucleates heterogeneously on the surfaces of C3S, so in this simulation the

nucleation barrier for C–S–H on C3S was set low enough to encourage C–S–H nucleation

exclusively on C3S surfaces. Therefore, a strong spatial correlation between C–S–H and C3S

is expected to arise due to nucleation effects. However, neither portlandite nor ettringite are

assumed to nucleate preferentially on any solid surfaces, so the locations of their precipitation

should be dictated only by the relative abundance of solute species from which they form.

Al(OH)

−

4 and H2SiO2−4 have lower intrinsic diffusion coefficients in water than Ca2+, OH−

or SO2−

4 , so one might expect to observe ettringite to form in the vicinity of C3A and C–S–H to

form in the neighborhood of C3S, the sources of Al(OH)

−

4 and H2SiO2−4 , respectively. In fact,

the preliminary results of this simulation indicate that the locations of ettringite and portlandite

are strongly correlated with C3A, but are not correlated significantly with either C3S or gypsum.

This can be seen visually in figure 5(c), where the proximity of ettringite crystals to tricalcium

aluminate is readily observed. To probe the effect quantitatively, radial distribution functions

(RDFs) for both ettringite and portlandite were calculated from the simulated microstructure

at 3 h of hydration. The RDF of a phase i gives the ratio of the concentration of that phase to its

average concentration, φi(r)/φi as a function of distance r from a reference point. Figure 7

plots the RDF for portlandite using surfaces of C3A, C3S or gypsum as the reference point.

For example, the curve with the circular points is the RDF using C3A surfaces as the reference

point, averaged over all of the C3A surface sites in the microstructure. Figure 8 gives the RDFs

calculated for ettringite using the same three surfaces as reference points.

Figures 7 and 8 show that both portlandite and ettringite are strongly correlated with

C3A surfaces. The concentrations of these phases near C3A are about 2.5 times their average

concentrations throughout the microstructure. Furthermore, the RDF of both phases decays

rapidly with distance, reaching a minimum value at about 8μm, where both phases are

12
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

2 4 6 8 10 12

r (μm)

0.0

0.5

1.0

1.5

2.0

2.5

3.0

φ / <φ>

C3A

C3S

Gypsum

Ca(OH)2

Figure 7. Calculated RDFs for portlandite using C3A, C3S or gyspum surfaces as the reference

point, for the microstructure shown in figure 5(b). The standard deviation of the sample mean at

each distance is less than 0.07, about the size of each point.

2 4 6 8 10 12

r (μm)

0.0

0.5

1.0

1.5

2.0

2.5

3.0

φ / <φ>

C3A

C3S

Gypsum

Ettringite

Figure 8. Calculated RDFs for ettringite using C3A, C3S or gyspum surfaces as the reference

point, for the microstructure shown in figure 5(b). The standard deviation of the sample mean at

each distance is less than 0.07, about the size of each point.

relatively depleted compared with their average concentrations. The location of the minimum

at 8μm corresponds closely to the average spacing of C3A domains in the microstructure, as

shown by plotting the autocorrelation function for C3A (figure 9). The first minimum of the

autocorrelation function is the average spacing of C3A domains [31]. Therefore, the increase

in the RDF at distances greater than 8μm is likely due to the approach of another C3A surface.

In contrast to their strong correlation with C3A surfaces, neither portlandite nor ettringite

are strongly correlated with C3S or gypsum surfaces; the RDF departs only slightly from unity

at any distance from these surfaces. The fact that ettringite is not strongly correlated with

gypsum indicates that SO2−

4 ions are sufficiently mobile that their diffusion does not limit

13
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

0 2 4 6 8 10 12

r (μm)

0.00

0.01

0.02

0.03

0.04

0.05

S2(r)

6 7 8 9 10

0.0027

0.0030

C3A

Figure 9. Calculated autocorrelation function for C3A for the microstructure in figure 5(b). The

inset shows a magnified view of the first minimum in the autocorrelation curve, which is interpreted

as the average spacing of C3A domains.

the formation of ettringite. Portlandite, which is a product of the hydration both of C3S and

of C3A, is not strongly correlated with C3S because silicate ions are not required to form it.

In addition, volume exclusion due to the presence of C–S–H near the C3S surfaces seems to

prevent portlandite from forming as abundantly near C3S as it does near C3A.

These kinds of spatial correlations among phases are now being investigated

experimentally using scanning electron microscopy coupled with energy dispersive

spectroscopy (EDS) to test the predictions made by HydratiCA [32]. A strong coupling of

experiment and simulation will help determine what, if any, general statements can be made

about the reaction network topology and spatial correlations in microstructures as complex as

those found in hydrating cement pastes.

4. Summary

A lattice-based cellular automaton model has been described for simulating 3D microstructure

development under multiple coupled nonlinear reactions and diffusive mass transport. The

algorithms applied at each lattice site require information only from the site itself and its

neighboring lattice sites. The localized nature of the calculations enables the model to be

parallelized to run on multiple processors. For the model system size investigated here, the

simulation speed scales well up to 100 processes, but the scaling behavior diverges significantly

from linear for higher process numbers. Likewise, the computational efficiency decreases

nearly linearly up to 448 processes, where inter-process communications reach 65% of the

total process time.

HydratiCA is a general model that should be applicable to a diverse range of reactive

systems, including porous geochemical environments and biochemical phenomena inside

living cells. Here we have focused on systems relevant to hydrating cementitious materials,

providing several points of verification and validation of the implementation on aqueous

suspensions of portlandite and ettringite. The model accurately captures the effects of

temperature and pH on the kinetics and equilibria of the dissociation of minerals such as

portlandite in water. In at least one instance, the simulations give significantly better agreement

14
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

with experimental measurements than the PHREEQC software package using the LLNL and

WATEQ4 databases. More chemically complex minerals such as ettringite, with exceedingly

small solubility products, are simulated with an accuracy that is comparable to the range of

measured or calculated solubilities that have been previously reported.

The coupling of realistic chemical kinetics and thermodynamics to 3D microstructure

enables HydratiCA to directly incorporate the effects of particle size distribution, particle

shape and surface area into its calculations. In addition, the model can be used to investigate

the emergence of spatial correlations between growing and dissolving phases in complex 3D

microstructures such as those found in hydrating cement pastes. Preliminary results for an

idealized Portland cement paste were shown in this paper, but more systematic experimental

and computational investigations are currently underway.

Minerals such as portlandite, gypsum and ettringite tend to have extremely anisotropic

growth rates in cement paste, leading to crystal morphologies that are tabular or acicular with

aspect ratios approaching 100. The current implementation of HydratiCA does not account for

growth anisotropy in any way, but it could do so by making the growth probability a function

of the unit normal vector to the interface at each lattice site where growth occurs. The form of

this function could be taken from polar plots of the so-called growth mobility function which

is used in the geometric theory of crystal growth [33].

Acknowledgments

This work was supported by the Virtual Cement and Concrete Testing Laboratory Consortium

and by the Partnership for High-Performance Concrete Technology program (HYPERCON)

at the National Institute of Standards and Technology. Large-scale parallel computations were

made possible by a grant from NASA/NAS. The author is indebted to Barbara Lothenbach

for insightful comments about thermodynamic calculations of solubility and to Ken Snyder,

Ed Garboczi and Jeff Thomas for thoughtful critiques of the manuscript.

References

[1] Bullard J W 2007 Approximate rate constants for nonideal diffusion and their application in a stochastic model

J. Phys. Chem. A 111 2084–92

[2] Bullard J W 2007 A three-dimensional microstructural model of reactions and transport in aqueous mineral

systems Modelling Simul. Mater. Sci. Eng. 15 711–38

[3] Bullard J W 2008 A determination of hydration mechanisms for tricalcium silicate using a kinetic cellular

automaton model J. Am. Ceram. Soc. 91 2088–97

[4] Bentz D P 1997 Three-dimensional computer simulation of cement hydration and microstructure development

J. Am. Ceram. Soc. 80 3–21

[5] Bishnoi S and Scrivener K L 2009 μic: a new platform for modelling the hydration of cements Cem. Concr.

Res. 39 266–74

[6] Karapiperis T and Blankleider B 1994 Cellular automaton model of reaction-transport processes Physica D

78 30–64

[7] Langmuir D 1997 Aqueous Environmental Geochemistry (London: Prentice-Hall)

[8] TaylorHFW 1997 Cement Chemistry 2nd edn (London: Thomas Telford)

[9] Parkhurst D L 1995 User’s guide to PHREEQC—a computer program for speciation reaction-path, advectivetransport, and geochemical calculations Water-Resources Investigations Report 95-4227 US Geological

Survey

[10] Mills R and Lobo V M M 1989 Self-Diffusion in Electrolyte Solutions (Amsterdam: Elsevier)

[11] Hummel W, Berner U, Curti E, Pearson F J and Thoenen T 2002 Nagra/PSI Chemical Thermodynamic Data

Base 01/01 (Parkland, FL: Universal Publishers)

[12] Lothenbach B and Winnefeld F 2006 Thermodynamic modelling of the hydration of Portland cement Cem.

Concr. Res. 36 209–26

15
Modelling Simul. Mater. Sci. Eng. 18 (2010) 025007 J W Bullard et al

[13] Pacheco P 1996 Parallel Programming with MPI (San Francisco, CA: Morgan Kaufmann)

[14] Tadros M E, Skalny J and Kalyoncu R S 1976 Kinetics of calcium hydroxide crystal growth from solution

J. Colloid Interface Sci. 55 20–4

[15] Lasaga A C 1981 Rate laws of chemical reactions Kinetics of Geochemical Processes (Reviews in Mineralogy

number 8) ed A C Lasaga and R J Kirkpatrick (Washington, DC: Mineralogical Society of America) pp 1–68

[16] Lambert I and Clever H L (ed) 1992 Alkaline Earth Hydroxides in Water and Aqueous Solutions (Solubility Data

Series) vol 52 (Oxford, UK: International Union of Pure and Applied Chemistry, Pergamon Press)

[17] Kulik D A 2002 GEMS-PSI 2.0 http://les.web.psi.ch/Software/GEMS-PSI/

[18] Thomas J J, Rothstein D, Jennings H M and Christensen B J 2003 Effect of hydration temperature on the

solubility behavior of Ca-, S-, Al-, and Si-bearing solid phases in Portland cement pastes Cem. Concr. Res.

33 2037–47

[19] Lea F M 1970 The Chemistry of Cement and Concrete (Washington, DC: Edward Arnold) 3rd edn

[20] Ghorab H Y and Kishar E A 1985 Studies on the stability of calcium sulfoaluminate hydrates: 1. Effect of

temperature on the stability of ettringite in pure water Cem. Concr. Res. 15 93–9

[21] Atkins M, Glasser F, Kindness A, Bennett D, Dawes A and Read D 1991 A thermodynamic model for blended

cements DoE Report DoE/HMIP/RR/92/005 US Department of Energy

[22] Zhang F, Zhou Z and Lou Z 1997 Solubility product and stability of ettringite Proc. 10th Int. Congress on the

Chemistry of Cement (Goteborg, Sweden)  ̈ pp II–88–I–93

[23] Alon U 2003 Biological networks: the tinkerer as an engineer Science 301 1866–7

[24] Alon U 2007 Simplicity in biology Nature 446 497

[25] Ridge M J 1964 Hydration of calcium sulphate hemihydrate Nature 204 70–1

[26] Gartner E M, Young J F, Damidot D A and Jawed I 2002 Hydration of Portland cement Structure and Performance

of Cements 2nd edn ed J Bensted and P Barnes (New York: Spon Press) pp 57–113

[27] Garrault-Gauffinet S and Nonat A 1999 Experimental investigation of calcium silicate hydrate (C–S–H)

nucleation J. Cryst. Growth 200 565–74

[28] Garrault S and Nonat A 2001 Hydrated layer formation on tricalcium and dicalcium silicate surfaces:

experimental study and numerical simulations Langmuir 17 8131–8

[29] Garrault S, Finot E, Lesniewska E and Nonat A 2005 Study of C–S–H growth on C3S surface during its early

hydration Mater. Struct. 38 435–42

[30] Thomas J J 2008 A new approach to modeling the nucleation and growth kinetics of tricalcium silicate hydration

J. Am. Ceram. Soc. 90 3282–8

[31] Berryman J G and Blair S C 1986 Use of digital image analysis to estimate fluid permeability of porous materials:

application of two-point correlation functions J. Appl. Phys. 60 1930–8

[32] Snyder K A 2009 Personal communication

[33] Taylor J E, Cahn J W and Handwerker C A 1992 Overview 1. Geometric models of crystal growth Acta Metall.

Mater. 40 1443–74

16