A Three-Phase Transport Model for High-Temperature Concrete Simulations Validated with X-ray CT Data

Concrete exposure to high temperatures induces thermo-hygral phenomena, causing water phase changes, buildup of pore pressure and vulnerability to spalling. In order to predict these phenomena under various conditions, a three-phase transport model is proposed. The model is validated on X-ray CT data up to 320 °C, showing good agreement of the temperature profiles and moisture changes. A dehydration description, traditionally derived from thermogravimetric analysis, was replaced by a formulation based on data from neutron radiography. In addition, treating porosity and dehydration evolution as independent processes, previous approaches do not fulfil the solid mass balance. As a consequence, a new formulation is proposed that introduces the porosity as an independent variable, ensuring the latter condition.


Introduction
Exposing concrete structures to fire can lead to spalling, a process where material violently breaks off of the surface. How susceptible a concrete mixture is to spalling is therefore an important criterion for the safety and service life of concrete structures. The prevailing hypothesis is that spalling is caused by a combination of high pore pressures and thermal stresses. When concrete is exposed to high temperatures, dehydration and free water evaporation often exceed the rate of vapour migration, resulting in rising pore pressures. Due to its low permeability, this is especially severe for high performance concrete (HPC), making it generally more susceptible to spalling. With a numerical analysis of the thermo-hygral transport in porous media in general, and of concrete in particular, insights into the mechanisms and influencing factors can be gained.
Many different approaches have been proposed [1][2][3], where the most accurate are the three-phase models [4][5][6], i.e., models that consider liquid water, moist air and the solid matrix as separate phases. These are usually validated using temperature and pressure gauges at various positions within the specimen. The pressure measurements in particular can be problematic. The presence of the sensor influences the resulting cement matrix around it, and localized cracks are more likely to develop in its vicinity [7]. Furthermore, the size of the sensor is much larger than typical pore diameters, resulting in an averaged value over the measured area. To overcome these experimental difficulties, non-destructive computer tomography measurements are increasingly used. Using these CT results for validation of the thermo-hygral transport model is an important improvement pursued in this work.
Balázs et al. have used X-ray CT to measure the change in density of the concrete specimen after undergoing a heating cycle [8]. A reduction in density by 2% to 4% was found after heating to 500°C, with good agreement between CT and conventional methods.
The use of CT can resolve these changes in space. In a separate fire test, the reduction in density primarily occurred close to the surface that was subjected to the fire. Henry et al. have used X-ray microtomography to characterize the change in total pore space and its connectivity after heating to 600°C for an hour [9]. The total pore space was found to be roughly 1.5 times larger, while the connectivity had increased by 3.5 times. The influence of air and water re-curing on the pore space, as well as the pore size distribution, was also investigated. In contrast to the experiments in [8,9], data from Powierza et al. [10] was obtained in-situ during the high temperature exposure. This makes it ideal for validating the water content in the pores to a numerical model and is, therefore, used here. The measured water content shows that the majority of previous models do not accurately reflect the water transport inside the concrete, as already pointed out in [7]. Additionally, knowing the water content allows for a deeper investigation of the influence of different dehydration formulations, several of which are compared in Sections 2.1 and 3.2. We are unaware of any such comparison in the literature. The suitability of thermogravimetric analysis as the basis for dehydration models will be discussed. There is a mass balance equation associated with every phase of the porous medium. Previous approaches did not solve the skeleton mass balance equation, assuming the error to be negligible. The novelty in this work is that the skeleton mass balance is solved along with the other balance equations, avoiding this error and resulting in one fewer constitutive equation. This also allows for dehydration descriptions that depend on the pore state, e.g., liquid water saturation or pore pressure, a significant extension compared to previous models. These modifications allow solving the coupled problem more accurately, while introducing the possibility of more complex dehydration models in the future.
The material is modelled as an isotropic homogeneous continuum. As a macroscopic model, it does not resolve the aggregates as discrete particles in the finite element mesh. The focus lies on the heat and moisture transport phenomena; coupling to mechanical deformation and damage evolution are possible extensions in future work. In the model, the aggregates are assumed to be inert, both chemically and physically, i.e., they release no dehydration water. Thermogravimetric analysis showed that basalt, granite and sandstone lose up to 3.1% of mass when heated to 1000°C [11], with similar results reported in [12]. Above 600°C, carbonate aggregates (limestone and dolomite) decompose into CaO and CO 2 [12], resulting in additional mass loss. Concretes with four different aggregates exposed up to 1000°C showed similar temperature-dependent tensile and compressive strength decrease [12]. This indicates that cement paste presents the weakest link in concrete. However, different aggregates contribute to microcrack evolution due to the thermal expansion mismatch and this will impact permeability, sorption isotherms and diffusivity.
We have made a conscious effort to make the reproduction of this work as convenient as possible. All the necessary code and data are openly available, a container with all the necessary software is provided, and running all the required steps is fully automated. See Section 4 for details.

Dehydration and Porosity Descriptions from the Literature
The importance of taking dehydration and changes in porosity into account when modelling the behaviour of concrete at high temperatures is widely recognized in the literature. We will present a number of dehydration and porosity descriptions.
Dehydration is the chemical process whereby chemically-bound water from the hardened cement paste is released into the pore space as temperatures increase. Heating above 105°C releases non-evaporable water (hydration water), which occupies part of the gel pores or is chemically bound in chemical phases. How much water is released is usually determined experimentally by thermogravimetric analysis (TGA). Different aggregate types exhibit different temperature stability. When using TGA data to model dehydration, the assumption that mass loss equates to water release only holds for aggregate types that are chemically and physically stable at high temperatures, such as quartzite or basalt.
Below 800 K, most aggregates are chemically and physically stable. Most models in the literature are calibrated using TGA measurement data.
For better comparison, all models have been brought to the same form, where c is the cement content (kg m −3 ), ν is the mass ratio of maximally released water to cement content, and Γ is the dehydration degree (Γ ∈ [0, 1]). A common assumption is that water loss below 105°C is water evaporation, with dehydration remaining at zero; the constant T d = 378.15 K will be used in these cases. In the formulation proposed by Pesavento, ν = f s m, where f s is the stoichiometric factor and m is the ageing degree [6]. Only the dehydration degree is assumed to be a temperature-dependent sigmoid function. Gawin et al. use the same product formulation as Pesavento, but propose a cubic polynomial for the dehydration degree [13].
The formulation by Tenchev et al. is equivalent to using a temperature-dependent ν, i.e., m dehyd (T) = cν(T). The proposed formulation is a step-wise linear function. Dwaikat and Kodur use a very similar expression [14].
The description by Dal Pont and Ehrlacher is one of the few approaches that takes the kinetics of the dehydration process into account [15]. The dehydration evolution is described by the differential equatioṅ where τ is the characteristic time of mass loss and m eq dehyd is the amount of water created at equilibrium. The value for the characteristic time is given in [15] as three hours.
All of the previous approaches are based on thermogravimetric analysis curves. TGA data are highly influenced by sample size, heating rate and vapour pressure [16]. The data used in these approaches were obtained on very small samples, low heating rates and in dry conditions. Furthermore, the results of TGA also include the CO 2 released from carbonates.
These conditions do not necessarily apply to specimens tested for water migration and spalling risk under high temperatures. As a way of correcting for these discrepancies while keeping the function-of-temperature approach, Dauti et al. perform an inverse analysis on measured neutron radiography data [7]. The result is a much steeper curve (see Figure 1). This result, however, cannot be directly used in other simulations, since it reflects the sample size, heating rate and pore vapour state of their experiments. The given graph also does not extend beyond 450°C.
Since the authors Dauti et al. have not provided an analytic description of their dehydration model, a logistic function was fitted to the given graph to obtain a mathematical expression as, For the curve presented in their paper, the best fit is obtained with a = 0.8219, k = 0.0876 K −1 and T 0 = 578.1 K.
Porosity is the ratio of pore space to the overall volume, n = (V w + V g )/V, where V w and V g are the volumes of the liquid and gas constituents in the averaging volume element V. In cementitious systems up to 105°C, the porosity is composed mainly of capillary pores (diameters approximately 10 nm-10 µm) [17], physically adsorbed water covering surfaces of gel particles, entrapped and entrained air and internal pores within aggregates. The evaporable water comprises free water residing in the capillary pores plus physically adsorbed water. The state of water and detailed description of a colloid C-S-H model with densities can be found elsewhere [18].  Different authors provided temperature-dependent porosity evolution, briefly summarized below and in Figure 2. Dauti et al. describe the porosity as a linear function [7], The porosity evolution according to Tenchev et al. is an empirical function [4], noting that factors such as dehydration, chemical decomposition of aggregates, thermal strains and microcracking contribute to its increase. Their actual description assumes the porosity to stay constant below 100°C, three times the initial porosity above 800°C and a cubic polynomial in the range between T ∈ [100°C, 800°C]: The coefficients of the polynomial are chosen such that the porosity is C 1 continuous. Dal Pont and Ehrlacher [15] give the change in porosity simply as n = n 0 + 0.72 × 10 −3 m dehyd . In the publications by Gawin et al., the porosity is described as a simple linear function, with initial porosities given between 0.06 and 0.087, and the constant A n in the range 1.63 × 10 −4 K −1 to 1.95 × 10 −4 K −1 . A much larger range of initial porosities (between 0.0512 and 0.13) are reported in later papers [19,20]. They also note that dehydration is only one factor driving the change in porosity, preferring the above empiric relation to a direct, dehydration-only description. Dwaikat and Kodur obtain the change in porosity from the sum of volume fractions [14] n = n 0 + V dehyd , where V dehyd is the volume of the dehydrated water and assumed to be V dehyd = m dehyd ρ w . However, this is in conflict with the skeleton mass balance, as discussed further in Section 2.2.

Balance Equations
The general form of the mass balance equations, namely of dry air, water vapour, liquid water and the solid matrix, is where m π is the mass per unit volume, J π is the mass flux and r π is the mass source, each with respect to phase π = {a, v, w, s}. The indices denote dry air, water vapour, liquid water and skeleton, respectively. For the solid phase, the mass balance is given by with the skeleton mass per unit volume m s = (1 − n)ρ s . This equation is often not solved in other models [4,13]. Dehydration generally causes an increase in porosity and further chemical changes in the solid. The idea of a constant skeletal density has been used as a C-S-H bulk density description in the whole isotherm range [18]. In such a case, the C-S-H sheet maintains a constant skeletal density of 2850 kg m −3 , and packing of the globules provides the correct C-S-H bulk density for different water contents [21]. Hydrated cement paste is mainly composed of C-S-H, together with portlandite, ettringite and other minor phases. We may extend the idea of a constant C-S-H skeletal density to the whole cement paste and concrete, resulting in a constant skeletal density ρ s . Such a description simplifies the effect of dehydration. Note that concrete bulk density provides no information on the skeletal density. This means that only two of the quantities in the set {n, ρ s ,ṁ dehydr } are independent. For the cases discussed here, the mechanical deformation, and therefore v s , is kept zero. The dry air mass conservation is described by where m a is the dry air mass per volume of a porous medium, and J a is the dry air mass flux. The mass per unit volume of dry air is given by where n is the porosity, S g and S w the gas and liquid water saturation, respectively, and ρ a is the density of dry air. Free water is present in two different phases-as water vapour and as liquid water. This results in two mass balance equations, These are commonly summed up to remove the evaporation termṁ vap , resulting in a balance of total water, The enthalpy balance can be written as where (ρC p ) eff is the effective heat capacity of the concrete, including water and moist air contributions, C pw is the specific heat capacity of liquid water, C pg is the specific heat capacity of moist air and λ eff the effective thermal conductivity of the concrete. The vapour mass source iṡ

Constitutive Equations
The dry air and water vapour flows are composed of an advective and a diffusive part, where v g is the gas velocity. For binary mixtures such as moist air, the rates of mass diffusion are of equal magnitude and opposite direction for constant total concentration. Therefore the gas diffusion coefficient is D = D AV = D VA . Additionally, it is assumed to be constant. The liquid water flow is only advective, where v w is the liquid water velocity. The advection is described by Darcy's law, The intrinsic permeability k increases as the temperature rises. This is due to both increasing porosity caused by dehydration and crack creation as a result of thermomechanical stresses. Here, it is given by where k 0 is the intrinsic permeability at the reference temperature T ref and A k is a material constant. This is a modification of the formulation proposed in [19], which specifies further terms for capturing gas pressure and damage dependence. The relative permeabilities depend on the current saturation with several models having been proposed. Here, we are using a description given by Beneš and Štefan [22], with ψ = 0.05 − 22.5n. A common choice for the sorption isotherms, also adopted here, is the model proposed by Baroghel-Bouny et al. [23]. Based on fitting experimental data for four different concrete mixes, the following formulation is obtained: where a and b are material parameters. The hysteresis between sorption and desorption branches is neglected. Solving for the saturation gives The capillary pressure p c can be calculated from the Kelvin equation, where p vs is the saturation vapour pressure. The effective heat capacity (ρC p ) eff is taken as a weighted sum of the specific heat capacities of the individual components, The specific heat capacity of the skeleton according to [24] is where The heat capacities of liquid water, water vapour and dry air can be found in Appendix A. The dehydration enthalpy ∆H dehydr is assumed to be constant [22]. The effective thermal conductivity λ eff is the product of the dry thermal conductivity λ dry and a factor for the increased thermal conductivity of pore water [5] where the dry thermal conductivity depends on the temperature

Boundary Conditions
The temperature boundary conditions are either Dirichlet conditions or convective-radiative Robin conditions where h, and σ are the convective heat transfer coefficient, the surface emissivity and the Stefan-Boltzmann constant, respectively. The air pressure is prescribed via a Dirichlet condition Lastly, the boundary condition for the vapour transport is also a Robin condition, where β c is the convective water vapour transfer coefficient.

Numerical Approximation
The set of balance Equations (10), (11), (15) and (16) is solved via the finite element method. When solving only the air mass balance, water mass balance and enthalpy balance, the set of independent variables is {p v , p a , T}. In the case where, additionally, the skeleton mass balance is solved for, the porosity n enters as a fourth independent variable. The choice of independent variable is not arbitrary. Gawin et al. argue that the capillary pressure p c can be used, even beyond the critical point of water, by reinterpreting the resulting values as a substitute for the product of water potential Ψ and liquid water density ρ w [19]. In our experience, the choice of {p g , p c , T} leads to poorer convergence.
Solving the additional skeleton mass balance results in a thermodynamically consistent solution, and in addition allows for one fewer constitutive equation. The characterization of the material may then be easier and faster, or measurements of the corresponding material parameter-porosity in this case-can be used as additional verification.
The domain is discretized with mixed elements, with a Lagrange element for each component. The vector of degrees of freedom of the unknowns can be written as The polynomial order is the same for all subelements, and linear elements were found to suffice. Increasing the order or choosing different orders for the individual components is trivial due to the implementation in FEniCS [25]. Since the test function w is from the same function space, it can also be split into its components w π .
To discretize Equations (10), (11), (15) and (16), they are transformed into their weak forms. The residual for the dry air mass balance is where J a is the dry air flux from Equation (18). The water mass balance reads with the vapour flux J v and the liquid water flux given in Equations (19) and (20), respectively. The skeleton mass balance residual is Lastly, the residual for the enthalpy balance, Since the dehydration is assumed to be an irreversible process, the dehydration degree is an internal variable that only depends on the maximum temperature at each material point, Discretization in time is achieved by applying the Rothe method, resulting in a nonlinear set of equations The time integration is performed using an Euler backward method with adaptive time stepping.

Skeleton Mass Density
Porosity and dehydration evolution are dependent processes, see Equation (10). Their evolution influences air, water and heat transfer, and their interdependence should not be neglected as in the previous models under discussion [4,6,7,[13][14][15]24]. As a result, the skeleton mass density is often treated inconsistently. It could be argued that the skeleton mass density changes in such a way as to fulfil the mass balance for two independently chosen descriptions of the porosity and dehydration. However, all models assume it to be constant when it enters into other equations, such as the formulation for the heat capacity.
Given the skeleton mass balance and a pair of dehydration and porosity models, one can solve the balance equation for a material point to obtain the skeleton mass density. To do so, the skeleton mass balance ∂m s ∂t =ṁ dehydr can be rearranged into an ODE for the skeleton mass density, This gives the theoretical evolution of that density such that the skeleton mass balance is maintained and can be seen for all the models under discussion in Figure 3.
The resulting values are higher than one would expect. Dry cement mixed with dry quartzitic aggregates in the same ratio as the concrete in Section 3.2 would result in a density of about 2740 kg m −3 . By hydration and curing, and then further dehydration by heating, no increase in skeletal density is possible. This demonstrates the problematic consequences of independently choosing porosity and dehydration descriptions.  Solving the skeleton mass balance directly links the dehydration and porosity evolution in our case. We assume the skeletal density to be constant. The additional balance equation allows for one fewer constitutive equation, thereby either saving on experimental effort or providing additional data for validation. If the dehydration description solely depends on temperature, the balance equation becomes an ODE that could be solved separately. Discretizing the skeleton mass balance along with the other balance equations allows for more complex dehydration descriptions, depending not only on temperature but also on, for example, liquid water content or pore pressure.

Validation on a Slowly Heated Cylinder
The data for validation have been obtained from an X-ray computer tomography (CT) scan of a concrete sample at various times during heating, published in [10]. A cylindrical specimen (∅40 ×100 mm) of high-strength concrete was surrounded by a glass-ceramic shell and then wrapped in aluminium-silicate wool to hinder the exchange of moisture and heat along the lateral boundary. The mix of the concrete can be seen in Table 1. The concrete has a compressive strength of 104 MPa [26]. The specimens were stored under water for 28 days after casting. At the time of measurement, they were at least 90 days old, having been stored in a climate chamber of 20°C and 65% relative humidity between water storage and heating. Heating was applied using an electric heating element at the top surface, with a heating rate of 10 K min −1 for the first 28 min up to a maximum temperature of 320°C. This temperature was then kept constant for another 130 min. The specimen was placed on a spinning table to take multiple images from different angles, and then reconstruct a 3D volume representation from the projections. Taking a single image took roughly one second, and a full scan was comprised of 650 projections, resulting in about 10 min per scan. As a consequence, the values at the given times are averaged values of 10 min around this point in time. The CT scans result in grayscale images where the luminance corresponds to the density of the material. The difference in brightness between empty and fully saturated pores was taken to represent 100% moisture change. The moisture change can thus be quantified by the change in brightness between the initial image before heating and subsequent images during heating. Deformation was corrected for using digital volume correlation. Since embedded temperature sensors would introduce artefacts into the X-ray images, separate experiments with embedded thermocouples were performed to find the temperature distribution. The top 5 mm of the images were removed, because they exhibit cone-beam artefacts rather than actual material behaviour. Further details of the experimental setup and image correction can be found in [10]. The images from the CT clearly show the change in moisture as the drying front advances into the material, see Figure 4. Averaging over the width of the diameter gives a one-dimensional moisture distribution, which we will use for comparison to the finite element model.  The initial conditions are given as a constant vapour pressure corresponding to a relative humidity of 65%, a dry air pressure of ambient air pressure (101.325 kPa) minus the vapour pressure, a temperature of 295 K and an initial porosity of 7.2455% as in [10]. Inhomogeneous initial moisture conditions due to drying during storage were initially considered, but found to be of little influence. The material parameters used in the simulation can be seen in Table 2.
Additionally, the temperature has been measured with four temperature gauges along the axial direction. A comparison of these results with the simulation can be seen in Figure 5. For comparison to the CT data, the change in moisture content is computed by where m max dehydr is the maximally dehydrated water mass when the dehydration degree Γ reaches one. The model shows good agreement, as illustrated in Figure 6. Particularly for the early part (15 and 30 min) and towards the end of the experiment (135 min), the location of the drying front of the simulation coincides with the experimental results. For the times in between (45, 60, 90 min), the model predicts a slightly more advanced front than can be observed in the experiments. The moisture accumulation behind the drying front, i.e., the moisture clog, is slightly overestimated at 30 min, and slightly underestimated for later times. The steep gradient of the moisture content is captured well. In the colder interior of the specimen, an increase in moisture content can be seen in the CT data. This broader accumulation of moisture is not present in the model results. There are two possible causes for this. Firstly, the dehydration description based on Equation (3) releases basically no water below 200°C (the temperature range of this region), which leads to an underestimation of the available water. Secondly, the analysis and postprocessing to get from CT brightness changes to water mass loss may introduce a systematic error as the experiment goes on.  The simulation allows for further insight into the state of water inside the pores, as well as the evolution of the solid phase. In Figure 7, the saturation reaches almost zero on the exposed surface, with a steep gradient, especially for the first 45 min. The moisture clog shows up again here. The maxima in saturation and gas pressure occur at the same time, yet not at the same locations. The maxima of the gas pressure are in front of the highest water accumulation, with values exceeding 30 MPa. For mechanical loads on the skeleton, the pore pressure is the relevant quantity, which is a purely postprocessed quantity given by

Figure 7.
Resulting saturation, gas pressure, liquid water mass, porosity, gas pressure, water pressure and pore pressure for each of the measurement times.
Due to the slow heating and comparatively low temperatures, the maximum pore pressure is about 5 MPa. The water pressure shows very large negative values in the dry regions of the specimen, but due to the low saturation there, the resulting influence on the pore pressure remains in the order of that of the gas pressure. As a result, the pore pressure remains mostly negative, since the capillary effect of the liquid water is larger than the increase in gas pressure. The porosity exhibits an almost step-like evolution due to the similarly steep dehydration formulation. After 60 min, it no longer changes since the temperature does not rise any further.
Of particular interest for this work is the choice of dehydration formulation. As mentioned in Section 2.1, dehydration conditions for TGA will differ significantly from those in a specimen for spalling tests. The models have all been recalibrated to the TGA data for this particular concrete. This was achieved using a nonlinear least squares method. The result of the calibration can be seen in Figure 8, and the parameters are shown in Appendix B. Note that while the concrete under investigation has quartzitic aggregates, the comparatively low temperatures of the experiment allow for experiment and modelling to be repeated for other aggregate types without difficulties. The simulation was repeated for each of the dehydration models that have been presented. The change in moisture mass shows significant differences between the models, see Figure 9. In particular, none of the models based on TGA exhibit enough dehydration during the early heating phase to match the experimental results. The dehydration description from Tenchev et al. leads to a loss of convergence beyond 21 min. The values during the first 21 min of heating are shown for completeness.

High-Temperature Benchmark Problem
While the CT data from the previous example enable direct validation of the moisture content inside the specimen, it is limited to temperatures of about 320°C. To show that the proposed model also works beyond the critical point of water, and to allow comparison to other models, a common benchmark problem is presented here.
It was proposed by Tenchev et al. [4] and also solved in [24]. A cross-section of a concrete column exposed to fire is simulated. Symmetry in the y-direction makes this a one-dimensional problem. The surface of the column is exposed to fire according to the ISO-834 fire curve for one hour. The outside atmosphere has a pressure of 0.1 MPa, and a relative humidity of 80%. The initial conditions for the concrete are a temperature of 20°C, a gas pressure of also 0.1 MPa and vapour pressure equal to the saturation pressure, i.e., a relative humidity of 100%. The resulting temperature and gas pressure distributions after 10, 30, and 60 min can be seen in Figure 10. For the temperature distribution, a very close match is obtained. Temperature information at 10 min. is not in the original paper. The gas pressure distribution appears to have a steeper decrease after the peak for our simulations.  In the paper that proposed this benchmark [4], an equidistant timestepping scheme with a timestep size of 2 s was used. The subsequent paper by Davie et al. [24] reduced it to 0.5 s. Both provided little justification for the specific choice of step size. Due to the highly nonlinear nature of the fire curve, an adaptive timestepping scheme can potentially save a significant amount of computing time. The criterion for increasing the time step is the number of Newton-Raphson iterations. If there are fewer than four iterations needed to solve the set of nonlinear equations, the new time step is set to 1.5 times the old time step. Consequently, much larger timesteps are used. Even with a conservative starting timestep of 5 s, the whole hour is integrated in just 87 timesteps, compared to 1800 or 7200 steps, respectively. The mean timestep was 41.44 s, and the largest timestep was 85.43 s. To compute the error introduced by larger timesteps, a reference simulation with equidistant timesteps of 0.01 s was performed. Equidistant timestepping with a timestep of two seconds as in [4] results in 0.0358, 0.0642 and 0.00202 as the L 2 norm of the relative errors for p v , p a and T. In comparison, the adaptive scheme with the much larger timesteps also leads to larger errors, with 0.0506, 0.0927 and 0.00306, again for p v , p a and T. For the largest error, that of p a , using the adaptive scheme leads to an increase in the error of about three percentage points, while computational effort is reduced by a factor of about 20.
Solving the skeleton mass balance leads to a consistency between the porosity and dehydration evolutions. If they are chosen independently, as was the case with all the previous approaches, the results violate the skeleton mass balance. The change in skeleton mass consists of two terms, the change due to porosity change ρ s (n 0 − n), and the change due to dehydration m dehyd . The incompatibility between the two terms can be seen in Figure 11.
The example shows that the model is able to cope with the transition beyond the critical point of water. Using an adaptive timestepping scheme can drastically reduce computation time. The introduction of the skeleton mass balance into the set of equations to be solved avoids incompatibility between dehydration and porosity descriptions.

Replication
This work, including all the tables and figures, is intended to be easily and fully repeatable. Reproducibility and repeatability are important parts of cumulative science. For a general introduction to reproducible computational research, see [27]. To allow repetition, finding mistakes or discrepancies between text and source code, as well as extend the analysis to new models or new data, an approach that is straightforward to use has been implemented in this work. The following section will explain its components.
Firstly, all of the input data, source code, parameters and scripts should be publicly available. In our case, all of the data are available at Zenodo, via the DOI 10.5281/zenodo.4452800 (accessed on 27 August 2021). This includes all the simulation code, TGA data, plotting instructions and LATEX files.
Secondly, in order to reproduce the results, not only is it necessary to have the source code and inputs but also the environment under which the code was executed, that is, the exact versions of all the programs and libraries used. We have chosen to provide a reproducible environment (or at least a reasonable approximation) in the form of a Docker container with all the required dependencies installed. Not only does this allow anyone to rerun the code in the same computational environment but also makes it very convenient. There is no need to manually install all the software that was used, which can be cumbersome and difficult. In addition, the Dockerfile serves as a readable description of all dependencies and the environmental setup. The container is available on Dockerhub at christophpohl/paper.
Lastly, and most importantly, all the data processing should be fully automated. This way, inefficient and error-prone manual steps are avoided, and a clear chain of reproducible steps is obtained. From running the simulation, via creating the graphs, to compiling the final PDF, all programs are invoked using a build automation tool. If reproducibility is the only criterion, this is not strictly necessary. A slightly simpler shell script that runs all required steps in sequence would suffice. However, changes in any of the files would require that everything is run again from scratch, even if only changes in the graphs or wording of small sections were made.
The use of a build automation tool avoids this by encoding the dependencies between the steps, their dependencies and their outputs as a directed acyclic graph (DAG). A visual representation of such a graph can be seen in Figure 12. The build tool will then assemble the DAG, and when asked to build the final PDF, traverse it to find the tasks that are not up to date and run them. A task is not up to date when either its output is missing or one of its dependencies has changed since the last time it was run. That way, changes in any of the predecessors will always trigger a reevaluation of all necessary tasks, and no running of unnecessary steps is taking place.

Conclusions and Outlook
Dehydration and its influence on the solid skeleton are key processes in understanding the moisture transport in concrete at high temperatures. Data from thermogravimetric analysis, which is performed on ground up concrete in artificial atmospheres, are not suitable for deriving a constitutive equation of dehydration for compact specimens under high heating rates. Comparing the predictions of existing models to X-ray CT data reveals drastic discrepancies, mainly related to their choice of dehydration description. Additionally, an inconsistency is introduced in these models by treating dehydration and porosity evolution independently. Solving the additional balance equation of skeleton mass yields the porosity as a result of the dehydration. That way, there is one less constitutive equation, and experimental data on the porosity may be used for additional validation.
The dehydration description proposed by Dauti et al. which is used here is based on an inverse analysis [7], and as a result is only valid within the bounds of the experimental setup. Extension to higher temperatures requires further experiments. Driving the research in this field is the need for higher spalling resistance, particularly for high performance concrete. Therefore, coupling the preceding model to damage mechanics in the search for better spalling prediction is a natural next step. Furthermore, mesoscale modelling of this problem might reveal interesting effects, as the interfacial transition zone has a higher porosity and connectivity [28], whereas the aggregates act as a barrier to moisture migration, increasing the tortuosity of the transport path. Lastly, obtaining a set of parameters for the constitutive equations is laborious, even for just one concrete. A publicly accessible repository of material data for concrete would greatly benefit the research, not only in this field.  The molar masses for air and water are M a = 28.971 × 10 −3 kg mol −1 and M w = 18.0153 × 10 −3 kg mol −1 . The density of the moist air is simply the sum of the densities of its components, ρ g = ρ a + ρ v .
The specific heat capacity of water is [24] C w p = and C a p = a 3 T 3 + a 2 T 2 + a 1 T + a 0 (A10) The evaporation enthalpy can be computed with the Watson formula [30], which for water reads ∆H vap = 2.672 × 10 5 J kg −1 T crit − T For better comparison, all models have been brought to the same form, where c is the cement content (kg m −3 ), ν is the ratio of maximally released water to cement content, and Γ is the dehydration degree (Γ ∈ [0, 1]). A common assumption is that the dehydration below 105°C remains zero; the constant T d = 378.15 K will be used in these cases. Pesavento gives the dehydration degree as a function of temperature [6] Γ = 0 for T ≤ T d , 1 2 1 + sin 1 2 π(1 − 2 exp(k(T − T d ))) for T > T d .
Dwaikat and Kodur use a step-wise linear function as well [14], for T > T 1 . (A16) Again, T 1 is not independent, with T 1 = T d + 1 r . As mentioned in Section 2.1, the approach of Dal Pont and Ehrlacher uses a differential equation [15]. The term for the equilibrium dehydration mass can be brought to the form of Equation (A12), and the dehydration degree can then be written as