# Thermodynamic modelling and simulation of hydrothermal sintering mechanism: understanding densification and porosity characteristic evolutions

S. Bordère, Lucas Villatte, Marie Anne Dourges, Graziella Goglio, Catherine Elissalde

► To cite this version:

S. Bordère, Lucas Villatte, Marie Anne Dourges, Graziella Goglio, Catherine Elissalde. Thermodynamic modelling and simulation of hydrothermal sintering mechanism: understanding densification and porosity characteristic evolutions. Acta Materialia, 2025, 301, pp.121544. 10.1016/j.actamat.2025.121544. hal-05278871

HAL Id: hal-05278871

https://hal.science/hal-05278871v1

Submitted on 23 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Thermodynamic modelling and simulation of hydrothermal sintering mechanism: Understanding densification and porosity characteristic evolutions

S. Bordère<sup>a,b,*</sup>, L. Villatte<sup>c</sup>, M.-A. Dourges<sup>d</sup>, G. Goglio<sup>c</sup>, C. Elissalde<sup>c</sup>

## ARTICLE INFO

### Keywords:

- Hydrothermal sintering

- Pore size distribution

- Thermodynamic modelling

- Numerical simulation

- Pressure solution mechanism

## ABSTRACT

Hydrothermal sintering of silica powders has never exceeded 75% densification, with reference to the literature. The study of phenomena likely to occur when densification stalls is addressed combining experimental and modelling approaches. Several samples were obtained by performing interrupted tests during a densification cycle of a submicron silica powder carried out at 300 °C for 60 min under a uniaxial pressure of 350 MPa. Despite the stagnation of densification, an evolution of pore size distribution at constant pore volume was evidenced using Mercury Intrusion Porosimetry (MIP). The thermodynamic modelling of the pore space aimed to establish a diagram determining the degree of liquid saturation of the porosity at 300 °C as a function of the relative density of the powder and the mass fraction of water initially mixed with the powder. Moreover, a numerical modelling of the mechano-chemical coupling of solute transfer at the solid/liquid interface, was specifically developed for this study. The thermodynamic diagram showed that the initial hydration of 7-10 wt.% of the silica powder ensures liquid water in the porosity throughout the thermal cycle, with significant wetting of the pore surface at 300 °C and 75% relative density. Simulation of the pressure-induced dissolution and precipitation processes in a 2D porous medium fully saturated with liquid and subjected to uniaxial pressure showed a local redistribution of matter along the pore surface in relation to the anisotropic stress field. Such non-densifying mechanism could be identified in light of these modelling results.

## 1. Introduction

A wide range of techniques have been recently developed in the field of ceramic sintering, both at high temperatures with extreme heating rates (>$10\,000\,^\circ\text{C/min}$), as well as at very low temperatures ($<500\,^\circ\text{C}$) involving chemical reactivity and exploring new driving forces to activate sintering mechanisms [1-4]. Lowering sintering temperatures and/or reducing the overall sintering time through fast heating is a timely challenge with environmental, economic and scientific impacts [5]. The development of non-conventional sintering methods involves different driving forces impacting transport mechanisms and thus microstructural evolutions. In this context, the emergence of very low temperature solvent-assisted sintering methods ($<500\,^\circ\text{C}$) has provided major contribution to the development of new materials and composites,

in particular for Cold Sintering Process [6-12]. Hydrothermal Sintering (HS) has also provided significant insights in the field of low temperature sintering [13-19]. According to the theoretical analysis of sintering at low temperature [20], hydrothermal sintering (HS) is close to uniaxial cold sintering in that an uniaxial pressure is applied to a compact in the presence of solvent. The main difference between the two processes is that one is carried out in an open system (CSP), while the other is carried out in a closed system (HS). This distinction gives rise to different evaporation characteristics, with the HS process favouring the retention of the solvent in its liquid phase.

Both HS and CSP are innovative processes inspired by geological densification processes in which compaction is mainly driven by intergranular pressure solution creep [21-26]. More specifically, the similarity between the HS or CSP and geological densification lies in this

<footnote><sup>a</sup> Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400, Talence, France</footnote>

<footnote><sup>b</sup> Arts et Metiers Institute of Technology, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400 Talence, France</footnote>

<footnote><sup>c</sup> Univ. Bordeaux, CNRS, Bordeaux INP, ICMCB, UMR 5026, F-33600 Pessac, France</footnote>

<footnote><sup>d</sup> Univ. Bordeaux, CNRS, Bordeaux INP, ISM, UMR 5255, F-33400, Talence, France</footnote>

<footnote>* Corresponding author at: Univ. Bordeaux, CNRS, Bordeaux INP, I2M, UMR 5295, F-33400, Talence, France.</footnote>

<footnote>E-mail address: sylvie.bordere@u-bordeaux.fr (S. Bordère).</footnote>
application of uniaxial pressure to the powder and the presence of an aqueous solution in which the dissolution/precipitation processes of the solid phase induced by stress gradients take place. When the dissolution occurs in the inter-granular hydrated area and precipitation at the neck surface, with a solute diffusion mechanism between these two reactions, densifying sintering occurs.

The presence of an inter-granular hydrated phase is more likely than a liquid film for the high pressures involved at grain contacts. This has been discussed by Biesuz et al. [20]. However, quantitative kinetic studies of densification carried out on ZnO using creeping law models have shown diffusion mechanisms through the contacts between the particles. These studies have also revealed that diffusion mechanisms at the grain contacts may be concomitant with contact sliding, which are both facilitated by dissolution/precipitation mechanisms. These studies were carried out on CSP [27,28] and HS [29].

If ZnO is particulary well suited for sintering at low temperature, densification of silica is limited. Silica was one of the most studied materials by hydro/solvothermal sintering. Pioneer's investigations were conducted by K. Yanagisawa et al. through hydrothermal hot pressing of silica powder synthesized by the Stöber process with different particle size [30–32]. Later, A. Ndayishimiye et al. carried out an in-depth study of the mechanisms involved in the densification of nanometer grain size silica powder mixed with water used as solvent [33]. Both the contribution of water and the impact of different experimental parameters on silica densification, were explored in detail. The authors have highlighted densification mechanisms based on mechano-chemical effects (i.e. pressure solution creep) in synergy with complex chemical ones such as polycondensation [33]. In line with this work, L. Villatte et al. have further investigated the sintering of 300 nm grain size silica by hydrothermal sintering [34]. They focused on the control of the surface chemistry of SiO₂ particles and on understanding their interactions with the solvent during the various steps of the process (preparation of the green pellet, sintering). The authors evidenced the fast desorption kinetics of the solvent by evaporation, which occurs after the addition of a known quantity of water to the dry powder during the mixing/grinding stage. A desorption threshold value has been identified. Its quantification and control during the initial stage are decisive to guarantee a reliable protocol and reproducible results.

It is not straightforward to draw precise conclusions on the optimum parameters for silica sintering conditions from the studies reported in literature, as different experimental conditions have been used (devices, particle size, surface chemistry, quantity and nature of solvent, pressure ranges). However, strong trends emerge: i) a liquid phase (solvent) is required to activate the dissolution processes; ii) the optimum sintering temperature when considering water as a solvent lies in the vicinity of 300 °C; iii) the maximum relative density obtained when using water as a solvent is quite limited, in the range 72–78 %, showing that water is not the most appropriate solvent to reach full density [35].

Beyond these trends, studies of silica sintered by hydrothermal sintering raise some important questions. One of the remaining issues is first to precisely control and quantify the amount of water not only in the initial powder-solvent mixture but also during the overall thermal cycle. The closed system conditions of HS make possible to limit the evaporation of the solvent in order to maintain a liquid phase in porosity that largely wets the surface of the pores or fills them entirely. Given these conditions, the second question is whether the pressure-solution mechanism can occur at the particle/liquid interface of the pores, becoming dominant compared to the grain boundary transport during the stagnation of the relative density. Drawing a parallel with solid-state sintering, for which the evaporation/condensation mechanism at the solid/gas interface with diffusion of species into the vapour phase induces a non-densifying process [36,37], this dissolution/precipitation mechanism located at the pore surface with diffusion through the liquid phase could also generate a non-densifying mechanism by HS.

To understand the low relative density obtained for silica densification and highlight the dominant mechanisms involved, it is therefore

essential to model the critical amount of liquid water at different steps of the process to determine whether this amount can guarantee potential optimum densification. In addition, pore-scale numerical modelling of dissolution/precipitation processes at the solid/liquid interface could, in combination with experimental results, make a valuable contribution to this understanding.

A preliminary section will present the experimental results of a case study, the hydrothermal sintering of silica with water. In particular, the evolution of the relative density of the powder will be correlated with the pore size distribution. Next, a thermodynamic model of the pore space will be developed, consisting of calculating the degree of saturation of the porosity with liquid water at the sintering temperature and during densification. This is the first milestone before moving on to numerical modelling of the mechano-chemical coupling of species transfer at the solid/liquid interface, which drives the pressure solution process. After validating the proposed model, the dissolution/precipitation process at the solid/liquid interface of pores of a partially dense ceramic will be simulated. The simulation results on the evolution of pore volume and pore morphology will be compared with the experimental data. The main objective behind this global approach is to demonstrate that modelling provides a significant contribution not only to a better understanding of the mechanisms involved in hydrothermal sintering, but also with a perspective of optimization of sintering conditions.

## 2. Case study: densification and pore characteristic investigation of silica nanoparticle powder sintered by hydrothermal sintering

### 2.1. Experimental protocol

In this experimental section, the aim is to evaluate and monitor the densification of commercial submicrometric grain size silica (NanoCym (Scottsdale, Arizona, USA) – 300 nm grain size), during hydrothermal sintering. To this end, our approach consists in characterizing various ceramics obtained from interrupted tests during a densification cycle carried out at 300 °C for 60 min under a uniaxial pressure of 350 MPa. The temperature and pressure conditions have been set on the basis of several previous studies and correspond to the optimum conditions when water is used as a solvent [34,35]. The first sample considered (sample 1 Fig. 1) is the green compact obtained by uniaxial pressure (UP) (60 MPa - 5 min) followed by cold isostatic pressing (CIP) (200 MPa - 10 min). Then all the others green samples were prepared as follows: the silica powder was mixed with the addition of 10 wt. % water and grinded. A fast desorption of the solvent can be observed, as a result the

<figure><img src="image_3.png" /><figcaption>Fig. 1. Thermal cycle of the different interrupted tests (gray curve) and time-dependence of the relative density: evaluation by mercury intrusion ● and geometric measurement ■.</figcaption></figure>
threshold value of water present in the mixing at the initial step is estimated between 7 and 10 wt. % [34]. The green compacts were then also prepared by UP and CIP. For the sintering cycle, a uniaxial pressure of 350 MPa was applied at room temperature and maintained constant throughout the cycle using a manual press [14,15]. Samples 2, 3 and 4 were obtained on the heating ramp at 170, 225 and 280 °C under 350 MPa, respectively. Samples 5, 6, 7 and 8 were obtained on the dwell temperature at 300 °C and 350 MPa at $t = 0, 10, 30$ and 60 min, respectively. For safety concern, the interrupted tests are carried out by pressure quenching and removal of the samples from the autoclave at 150 °C. All sample preparation conditions are listed in Table 1. The reproducibility has been checked performing twice samples for each interrupted test.

## 2.2. Characterisation of apparent density and pore size distribution during densification

The sample experimental apparent density and pore size distribution were investigated by Mercury Intrusion Porosimetry (MIP), using an Autopore IV 9500 porosimeter (Micromeritics) with the following parameters: contact angle = 130°, Hg surface tension = 485 mN.m⁻¹. Pressures are applied from 0.0034 up to 160 MPa (pore size diameter 367 µm/8 nm). The apparent density ($ρ_{app}$) using MIP technique at the low pressure was compared with that obtained using geometrical measurements. These values were used to calculate the relative density (Fig. 1). The value of 2.14 g/cm³ was used as solid density ($ρ_{skel}$), corresponding to the theoretical density value of silica. The relative density (%) was evaluated for each sample considering the formula: $d_{rel} = ρ_{app}/ρ_{skel}x100$. The results from the interrupted tests are reported in Fig. 1. A significant increase of the relative density is observed as the temperature rises, reaching an average value of 73 ± 2 % at the beginning of the dwell. This value remains low and relatively constant throughout the dwell at 300 °C.

The porosity of the samples is found to be distributed into mesopores (between 8 and 50 nm in these experimental conditions) and macro-pores (>50 nm) as it is reported on Fig. 2 with the evolution of the cumulative intrusion volume obtained on the samples 5, 6, 7 and 8. From this cumulative intrusion volumes, the proportion of mesopores is determined for each sample. It corresponds to the ratio of the cumulative volume of pores between 8 nm and 50 nm in size to the total cumulative volume. This evolution of the mesoporous contribution to porosity is shown in Fig. 3. It is worth noting that while the relative density remains roughly constant during the dwell temperature (Fig. 2), the porosity size distribution evolved slightly but significantly towards an increase in the proportion of mesopores (Fig. 3). There is thus a change in microstructure without densifying mechanism.

## 3. Thermodynamic modelling of pore space characteristics

It should be remembered that all studies performed on the hydrothermal sintering of silica agree on two points: i) the presence of solvent in liquid form is necessary to activate the densification mechanisms, and

<caption>**Table 1**</caption>

Conditions for obtaining silica ceramics samples by the interrupted test method. All the green samples hydrated or not were prepared by uniaxial pressure (60 MPa - 5 min) followed by cold isostatic pressing (200 MPa - 10 min).

<table><thead><tr><th>samples</th><th>Hydrated green compact</th><th>Temperature (°C)</th><th>Dwell time (min)</th></tr></thead><tbody><tr><td>1</td><td>no</td><td>20</td><td>-</td></tr><tr><td>2</td><td>yes</td><td>170</td><td>-</td></tr><tr><td>3</td><td>yes</td><td>225</td><td>-</td></tr><tr><td>4</td><td>yes</td><td>280</td><td>-</td></tr><tr><td>5</td><td>yes</td><td>300</td><td>0</td></tr><tr><td>6</td><td>yes</td><td>300</td><td>10</td></tr><tr><td>7</td><td>yes</td><td>300</td><td>30</td></tr><tr><td>8</td><td>yes</td><td>300</td><td>60</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Cumulative intrusion volumes obtained using Mercury Intrusion Porosimetry, for silica samples sintered at 300 °C at different sintering times: Ⓐ 0 min; Ⓑ 10 min; Ⓒ 30 min; Ⓓ 60 min. Note that mesopores with diameters less than 8 nm are not considering with this method.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Mesoporous proportion (±2 %) of porosity determined by Mercury Intrusion Porosimetry (pore size in the range of 8 50 nm) for silica samples sintered at 300 °C at different sintering times: Ⓐ 0 min; Ⓑ 10 min; Ⓒ 30 min; Ⓓ 60 min.</figcaption></figure>

ii) the optimum densification temperature is of the order of 300 °C [34, 35]. The mass fractions of water considered for the present experimental study (of the order of 7–10 wt. %) lead to a densification of the powder at a relative density of about 73 %, thus suggesting the presence of liquid water. A general methodology based on the thermodynamic properties of water is proposed here. This methodology allows to establish the link between the mass fractions of water evaluated during the sample preparation and the saturation degree of liquid water in the pore space at the sintering temperature and during densification.

Regarding the pore space characteristics to be modelled and quantified, the different scenarios likely to occur during the hydrothermal sintering process are shown on Fig. 4. The pore volume contains a fluid domain which may be single-phase (only vapor or liquid phase) or two-phase (liquid and vapor).

The name of the variables and constants of the problem are defined in Appendix A, and the complete methodology for the pore space characteristics is described in Appendix B.

The liquid solvent saturation degree in porosity, which is the property of interest to evaluate, corresponds to the volume fraction of the liquid occupying the porosity volume. All liquid saturation degrees in porosity $f_L^T(d)$ at sintering temperature $T$ and powder relative densities $d$ will depend on: i) the initial liquid saturation $f_L^{T0}(d_0)$ at room temperature.
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Illustration of the different scenarios likely to occur during the sintering process with regard to the evolving characteristics of the fluid domain contained in the porosity. The first step shows the changes with increasing temperature from $T_0$ (room temperature) to $T$ at constant initial density $d_0$ leading to a single-phase or two-phase pore space. The second step at constant temperature $T$ and increasing relative density shows that the single-phase pore space can remain single-phase or be transformed into a two-phase pore space. The two-phase pore space remains two-phase until a critical relative density is reached where the pores are completely filled with liquid. Note that the dark blue corresponds to liquid water and the light blue to water in vapor phase. Two characteristic volumes are indicated, that of the cell containing the sample $V_C$, and that of the four reservoirs $V_R$. The function of these reservoirs is to recover the solvent that may be expelled during powder densification, when the pores are completely filled with liquid. Pore space and reservoir volume are connected.</figcaption></figure>

ture $T_0$ and initial relative density $d_0$; ii) the pure water phase diagram (see Fig. B.1); iii) the sintering temperature $T$, iv) the powder relative density $d$; device characteristics in terms of sample cell volume $V_C$, and reservoir volume $V_R$ (Fig. 4). Note that the initial liquid saturation $f_L^{T_0}(d_0)$ is deduced from the weight percentage of solvent $f_{mL}^{T_0}$ in the initial compact, and by considering the initial relative density of the powder $d_0$, as well as the densities of liquid water $\rho_L^{T_0}$ and silica $\rho_S^{T_0}$ at $T_0$ (see Eq. (B.2)).

Two assumptions are made. First, for each configuration, the entire fluid domain is assumed to be in thermodynamic and mechanical equilibrium, and liquid water, when present, is considered to be located at the surface of the particle. Secondly, for the configuration at room temperature $T_0$ and initial relative density $d_0$, the amount of water in its vapor form has been neglected in the calculation of the volume fraction $f_L^{T_0}(d_0)$ of liquid water in the porosity (Eq. (B.2)). This is justified in the case of water as a solvent, since the saturation vapour pressure is 2 % of atmospheric pressure (see Table 2 for physical constants of water). So, in the case of water, the pore phase in this initial state can be considered to consist mainly of air and liquid solvent.

<caption>**Table 2**</caption>

Values of saturation vapor pressure $p_s^T$, liquid density $\rho_L^T$ and vapor density $\rho_L^T$ at the vapor saturation point for different characteristic temperatures. These data were obtained from the REFPROP database [38].

<table><thead><tr><th>T (°C)</th><th>p<sub>s</sub><sup>T</sup> (MPa)</th><th>ρ<sub>L</sub><sup>T</sup> (kg·m<sup>-3</sup>)</th><th>ρ<sub>V</sub><sup>T</sup> (kg·m<sup>-3</sup>)</th></tr></thead><tbody><tr><td>20</td><td>0.0023177</td><td>998.19</td><td>0017.162</td></tr><tr><td>300</td><td>0.8588</td><td>712.14</td><td>46.17</td></tr><tr><td>325</td><td>12.051</td><td>654.33</td><td>70.506</td></tr><tr><td>350</td><td>16.529</td><td>574.71</td><td>113.61</td></tr></tbody></table>

### 3.1. Degree of liquid saturation in porosity at sintering temperature T and initial relative density $d_0$

As a first step, the change in fluid domain characteristics is modelled during the initial temperature increase from $T_0$ to the sintering temperature $T$ at constant powder relative density $d_0 = 52 \%$ (first step shown in Fig. 4). As discussed before, pore space is constituted of liquid water and air at ambient temperature $T_0$. So, from various liquid water saturation degree $f_L^{T_0}(d_0)$ of the initial porous medium, the liquid water saturation degree $f_L^T(d_0)$ was calculated for three sintering temperature values (300 °C, 325 °C and 350 °C) (see Appendix B). To sum up, the successive calculation steps and assumptions are:

i) calculation of the apparent fluid density $\rho_M^T(d_0)$ (Eq. (B.7)) corresponding to the water mass (whatever in liquid or gaseous form) distributed into the entire available volume including pore space and reservoir volume; this apparent density depends on the cell volume $V_C(d_0) = 0,4 \text{ cm}^3$ and the total reservoir volume $V_R = 0.5 \text{ cm}^3$ of the hydrothermal sintering device (see Fig. 4) through $V_L^{T_0}(d_0)$ (Eq. (B.4) and $V_F^T(d_0)$ (Eq. (B.6));

ii) comparison of the calculated apparent fluid density $\rho_M^T(d_0)$ with those of the liquid phase $\rho_L^T$ and vapor phase $\rho_V^T$ at the saturation point at $T$ that are determined from the REFPROP database [38] (Fig. B.1 and Table 2);

iii) if $\rho_V^T < \rho_M^T(d_0) < \rho_L^T$, the fluid domain is diphasic and characterized by the saturation vapor pressure $p_S^T$, which does not depend on the value of $\rho_M^T(d_0)$. Note, that for temperature range 300 °C-350 °C, the saturation vapor pressure is about two order of magnitude higher than the atmospheric pressure, indicating that vapor phase is almost entirely constituted of water. Liquid water saturation $x_L^T(d_0)$ relative to the total available volume (pore volume and reservoir volume) can then be determined from $\rho_M^T(d_0)$, $\rho_V^T$ and $\rho_L^T$ (Eq. (B.8)). The liquid water saturation degree of porosity $f_L^T(d_0)$ is deduced from $x_L^T(d_0)$ by
considering the pore volume to the total fluid volume ratio (Eq. (B.9));

iv) if $\rho_M^T(d_0) < \rho_V^T$, the fluid domain is a vapor single-phase ($x_L^T(d_0) = 0, f_L^T(d_0) = 0$) whose pressure depends on the value of $\rho_M^T(d_0)$.

The calculations of the thermodynamic properties of the pore space at T = 300, 325, and 350 °C are shown in Fig. 5. First, the variations of the pore space pressure as a function of the degree of the initial liquid water saturation in porosity $f_L^T(d_0)$ (full lines) clearly exhibit plateaus for the three sintering temperatures. These plateaus reflect two-phase domains, within the volume combining porosity and reservoir volumes, at thermodynamic saturation vapor pressure $p_T^s$. For pressure values lower than $p_T^s$, the pore space is single-phase, thus only constituted of water vapor. For instance, for $f_L^T(d_0) = 0.3$, the porosity is characterized by a two-phase domain for T = 300 °C and T = 325 °C, whereas it is characterized by a single-phase domain at T = 350 °C. Secondly, the value of the new liquid saturation degree in porosity at each sintering temperature can be determined from the dashed lines in Fig. 5. As expected, these saturations decrease as the sintering temperature increases. For examples, for $f_L^T(d_0) = 0.3$ the saturations are $f_L^{300}(d_0) = 0.2$; $f_L^{325}(d_0) = 0.078$ and $f_L^{350}(d_0) = 0$ for T = 300 °C, T = 325 °C, and T = 350 °C, respectively. While temperature generally contributes to increase sintering kinetics, it also contributes to evaporate all the liquid water contained in the pores, thus preventing from solid dissolution and consequently limiting also the hydrothermal sintering process.

## 3.2. Liquid saturation degree in porosity at sintering temperature T and density d

As a second step, the change in fluid domain characteristics during powder densification at constant temperature T from initial relative density $d_0 = 52 \%$ to dense material, is modelled. More precisely, for a specific temperature T, the evolution of the liquid saturation degree of the pores $f_L^T(d)$ is plotted as a function of the relative density d for different values of the initial saturation $f_L^{T0}(d_0)$. Fig. 6 shows changes in liquid saturation of the porosity for the sintering temperature of 300 °C,

<figure><img src="image_7.png" /><figcaption>Fig. 5. Diagram showing the thermodynamic state of water in porosity at different sintering temperatures (T = 300, 325, 350 °C) as a function of the degree of the initial liquid water saturation $f_L^{T0}(d_0)$ at room temperature $T_0 = 20$°C and initial relative density $d_0 = 52 \%$. As a function of $f_L^{T0}(d_0)$ is drawn in full lines the variation of the pore space pressure p (left axis), and in dashed lines the liquid saturation degree $f_L^T(d_0)$ at T and $d_0 = 52 \%$ (right axis).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Evolution of the degree of liquid water saturation in porosity at T = 300 °C as a function of relative density d and for different initial saturations at room temperature $f_L^{T0}(d_0)$ (%): 50; 40, 30, 20, 15, 12.5. The grey area represents the saturation range in which the initial hydrated silica powders (sample 1 with 7–10 wt. % of water) are expected to evolve during densification.</figcaption></figure>

that can be discussed as follows:

i) given the assumption that the liquid phase is preferentially located at the surface of the particles (no liquid phase is present in the reservoirs until the porosity is fully saturated), total saturation of the pores can be achieved during densification, with the two-phase domain occupying only the volume of the reservoirs in this case.

ii) for each initial saturation degree $f_L^{T0}(d_0)$ considered, the liquid saturation degree values obtained at density $d_0 = 0.52$ and temperature T = 300 °C correspond to that already plotted in Fig. 5. For example, for liquid saturation degree of 30 %, the value $f_L^{300}(d_0) = 0.20$ is obtained.

iii) as shown in Fig. 4 and Fig. 5, the single-phase pore space can remain single-phase. The calculations show that this arises for initial saturation values lower than 12 %. The single-phase pore space can also be transformed into a two-phase pore space. This is the case for initial saturation values of 12.5 % and 15 % (Fig. 6), the formation of the liquid phase appearing for relative density around $d = 0.9$ and $d = 0.6$, respectively. Once the liquid phase is present, the saturation value increases during densification until the liquid phase occupies the entire pore volume $f_L^{300}(d) = 1$. Full saturation of the porosity may occur just before full densification of the powder, as in the case of initial saturation $f_L^{T0}(d_0) = 12.5 \%$, or may occur at low relative densities, such as $d = 0.67$ for $f_L^{T0}(d_0) = 50 \%$.

Given the initial relative density of the powder $d_0 = 0.52$, and the mass fraction of water introduced, which is estimated to be in the range $0.07 < f_{mL}^{T0}(d_0) < 0.1$ (see Section 2.1), the initial saturation of liquid water in porosity $f_L^{T0}(d_0) \times 100$ is assessed to be within the range [18–27 %] (Eq. (B.2)). This liquid saturation can be roughly converted into a ratio $e/r$ between the thickness of the water layer, e, and the radius of the particles r. It is defined approximately before densification and for small ratio by (see Appendix C):

$$ \frac{e}{r} = \frac{f_L^T(d_0)(1 - d_0)}{3d_0} \qquad (1) $$
So, at room temperature and initial relative density $d_0 = 0.52$, the ratio $e/r$ lies in the range: $0.055 < \frac{e}{r} < 0.083$. It gives rise, for particle radius $r = 150$ nm (grain size of experimental silica powder is 300 nm), to a water layer thickness in the range [8 – 12 nm]. By considering as 2.82 Å the mean van der Waals diameter of water [39], the number of water molecule layers around particles is roughly in the range [28 – 42 layers]. Note that for lower submicrometric grain size (100 nm), the powder specific area increases leading to a number of water molecule layers third times lower. Therefore, reducing particle size to increase sintering efficiency can, for a given mass fraction of water, lead to the opposite result.

During the temperature increase until $T = 300$ °C, the water evaporates leading to a liquid saturation range $0.021 < f_L^T(d_0) < 0.15$ in the pore space for $d_0 = 0.52$ (see Fig. 6). This liquid saturation decrease is counterbalanced by the powder densification during the heating ramp to reach the saturation range $0.18 < f_L^T(d) < 0.50$ for $d = 0.75$ (Fig. 6). In other words, on the basis of the thermodynamic calculations, the mean rate of evaporation of the liquid water initially present in the porosity (7 – 10 wt. %) can be estimated at 65 % when the temperature rises to 300 °C at a relative density of 52 %. For densification from 52 % to 75 % relative density, around 27 % of this evaporated water re-condenses. Then, about 50 % of the water initially introduced into the powder would be maintained in liquid form during the non-densifying stage.

This means that for the experimental conditions and densification results described in Section 2, the pore space of samples 5 to 8 (Fig. 1) is not fully saturated with water, but remains two-phase at the saturation vapor pressure. Nevertheless, it can be expected that the water layer thickness is large enough for the silica solute to diffuse through it from one point located on the solid/liquid surface to another.

## 4. Numerical modelling of dissolution/precipitation processes involving solid/liquid interface

The 3D direct simulations of the overall HS process (volume shrinkage, pore coalescence and closure of interconnected porosity) imply the resolution of a large number of physical equations describing several mechanisms: i) species transport at the grain boundaries and in the fluid, ii) pressure-induced dissolution and precipitation processes at the interfaces, iii) the mechanical liquid/solid interaction related to interface displacements, iv) the liquid and vapor flows in porosity even induced by capillary effects, v) the mass transfer at the liquid/vapor interface when liquid evaporates during temperature increase and vapor condensates during pore volume shrinkage. In the framework of this study, such 3D simulations are currently unrealistic.

This is why we have focused on exploring the apparently contradictory observation that, although not densifying, sintering is nevertheless associated with an evolution in the pore size distribution of open porosity on the dwell at 300 °C, as if densification was continuing. This suggests that the pore surface contributes to mass transfer through dissolution and precipitation phenomena. We therefore propose an original numerical model dedicated to this specific problem of pressure solution mechanism at the solid/liquid interface of pores. We expect this simulation approach to help understand the evolution of pore size distribution at constant pore volume.

Note that from the previous assessment of liquid water saturation at $T = 300$ °C and $d = 0.75$, the liquid layer along the solid/liquid surface was shown to be sufficiently thick for the silica solute to diffuse through. This is the necessary condition to activate the pressure-induced dissolution and precipitation processes at the pore surface.

In this numerical modelling section, we will first describe the physics of the mechano-chemical coupling at the interface which drives the pressure solution phenomenon (Section 4.1). Then the equations relative to the proposed modelling methodology of mass transport and transfer at the interface will be presented in Section 4.2. After methodology validation (Section 4.3.1), the simulation of a 2D porous media

uniaxially pressed will be presented and analysed in Section 4.3.2.

### 4.1. Phenomenological description of pressure solution phenomenon at solid/liquid interface

The pressure solution phenomenon is described on the basis of a two-phase system, for which the solid and the solution are referred to $\beta$ and $\alpha$ phases, respectively. Here, we consider that the $\alpha$-phase is a binary A-B solution, where A is the solvent (water) and B the silica-solute. The interface between the solution and the solid is noted $\Sigma$. Note, that all the complex chemistry of the dissolution of silica in water is here neglected. The $\beta$-phase is a pure phase (silica) only constituted of B-species. Let's consider the saturation concentration $c_B^{a*}$ of the B-solute in the solution, at reference pressure $p^*$ and temperature $T$, with $p^*$ identical for the liquid and the solid. If the pressure at the interface on the side of the solid increases, the induced pressure jump $(p_\Sigma^\beta - p^*)$ at the interface solid/solution, leads to the shift of the solute concentration at the interface in order to maintain the mechano-chemical potential constant through the interface. The equilibrium concentration $c_{B\Sigma}^a$ of the B-solute at the interface resulting from the mechano-chemical coupling is defined as a function of $(p_\Sigma^\beta - p^*)$ and $c_B^{a*}$ through the relation [40–42]:

$$c_{B\Sigma}^{a} = c_{B}^{a*} e^{\frac{\Omega(p_{\Sigma}^{\beta} - p^{*})}{RT}}, \quad (2)$$

where $\Omega$ is the molar volume of the B-species. It is considered identical for $\alpha$-phase and $\beta$-phase and independent of pressure, which is an acceptable assumption for condensed liquid or solid phases even with high pressure values.

In Fig. 7 is shown a schematic representation of the involved mechanisms during hydrothermal sintering. The pressure applied to the powder compact leads to a stress field within the volume of the grains, the pressure of the liquid for a two-phase pore space corresponds to the saturation pressure $p_s(T)$ at temperature $T$ (see Appendix B). In this model, the reference pressure $p^*$ corresponds to this saturation pressure and remains constant during densification (see Section 3.1). These volumetric stress fields emerge at the surface of the grains, creating a surface stress gradient. Stresses generated by surface tension are considered negligible compared to those generated by applied loading

Thus, the pressure jump $(p_\Sigma^\beta - p^*)$ varies along the interface inducing solute concentration gradient $\nabla c_{B\Sigma}^a$ along this interface, as well. Consequently, diffusion transport of the B-species through the solution is activated. This is illustrated in Fig. 7 where some points at the interface characterized by pressure $p_{\Sigma1}^\beta$ and concentration $c_{B\Sigma1}^a$ then have an undersaturated solution in their vicinity, leading to dissolution of the solid. Conversely, other points characterized by pressure $p_{\Sigma2}^\beta$ and concentration $c_{B\Sigma2}^a$ have a supersaturated solution in their vicinity, leading to precipitation of the solid (Fig. 7). To sum up, the pressure gradient at the interface induces a sequential mechanism of solid dissolution, solute diffusion and solid precipitation. It should be noted that this sequential mechanism may be located in the liquid film in-between two grains giving rise to densifying sintering [22,26,43], or at the liquid/solid interface in the pore space.

### 4.2. Mechano-chemical modelling of mass transport and transfer at the interface

The dissolution/precipitation phenomenon is complex to model and simulate at the interface scale since it requires to deal with concentration jumps at the interface, as well as moving interfaces. This applies to the scientific communities in the fields of computational fluid dynamics or computational materials science. In computational fluid dynamics, the main challenge is to deal with the phenomena of bubble dissolution and precipitation in liquid flows where, in addition to jumps in species
<figure><img src="image_9.png" /><figcaption>Fig. 7. Schematic representation of the sequential dissolution/diffusion/precipitation processes involved in hydrothermal sintering. In the solid, the pressure gradient resulting from the uniaxial pressing is illustrated by two pressure value $p_{\Sigma1}^{\beta}$ and $p_{\Sigma2}^{\beta}$ with $p_{\Sigma1}^{\beta} > p_{\Sigma2}^{\beta}$. The solution is characterised by the saturation pressure $p^{*} = p_{s}(T)$ of the two-phase pore space (see Section 3.1 and Appendix B). At the interface, pressure jumps between the solution and the solid ($p_{\Sigma1}^{\beta} - p^{*}$, $i = 1, 2$) displace the saturation concentration of the solute from $c_{B\Sigma1}^{\alpha}$ to $c_{B\Sigma2}^{\alpha}$, as formalized by Eq. (2). The larger the jump, the greater is this shift. This leads to impose conditions of $i$) undersaturation at the point where the pressure jump is high as illustrated by the concentration profile ①, and ii) oversaturation where the jump is low as in profile ③. Such profiles lead to dissolution ① and precipitation processes ③, respectively and transport of the solute in the solution ② from interface point ① towards interface point ③.</figcaption></figure>

concentration, there are also significant jumps in fluid density at interfaces [44–48]. In computational materials science, the key point is to model such dissolution/precipitation phenomena in solid/solid elastic media for morphology evolution of precipitates, where pressure jump at the interface is coupled with chemical equilibrium [49–54]. For solid/liquid systems with moving interfaces, it is also the interaction between the fluid and the solid having different rheological behaviours which is complex to deal with. Such solid/liquid interaction involving interface reaction is modelled in porous media problems [55–57]. The simulation of hydrothermal sintering combines all these numerical issues since it requires to deal with concentration and pressure jumps at the interface, as well as moving solid/liquid interfaces. Simulations of Ostwald ripening induced by surface tensions address these different points [58,59] and have made it possible to reproduce the coarsening of particles in the final stage of liquid-phase sintering. A recent article covers this subject to reproduce the pressure solution creep at the intergranular region induced by applied pressure [60].

### 4.2.1. Methodology

The proposed model for species transport and their transfer at the interface under non-isobaric condition is based on a one-field approach developed in computational fluids mechanics, meaning that a single equation describes the physics within the two phases and at the interface [47–48,57]. In comparison, in two-field approaches, each phase is characterized by its own set of equations with additional equations at the interfaces to respect mass conservation [44–46]. This modelling also differs from phase field modelling largely developed in computational material science, and which is based on a free energy functional [51–53, 55,58–60].

The mechanical modelling of the elastic stresses both in the solid induced by the pressure applied to the compact, and of the liquid flow

induced by the moving interfaces, is also described using a single-field approach. The mechanical equation of momentum conservation has already been developed for solid/solid and solid/liquid systems [57,58]. It has been shown in the latter case that it unifies under a single formulation the isotropic elastic behaviour of a solid with the viscous compressive flow of a fluid [61,62].

In one-field approaches, a phase function $\phi$ is considered to differentiate between the two phases with $\phi = 0$ for the $\alpha$-phase (solution) and $\phi = 1$ for the $\beta$-phase (solid). This phase function allows to define any physical constants $Y$ of the problem at every points of the system from the constant values $Y^{\alpha}$ and $Y^{\beta}$ of each phase through the relation:

$$Y = (1 - \phi)Y^{\alpha} + \phi Y^{\beta}. \qquad (3)$$

$Y$ is defining either the $k$-species diffusion coefficient $D_k$, the saturation concentration of the $k$-species $c_k^*$ at reference pressure $p^*$ and temperature $T$, the compressibility coefficient $\chi_T$ or the Lamé's coefficient $\mu_E$. The solid/liquid-interface is represented by the value of the phase function $\phi = 0.5$.

### 4.2.2. Mass transport and transfer equation

The single equation for diffusion species transport and transfer at interfaces was first developed for isobaric conditions [63] and then extended to non-isobaric conditions in the context of phase transformations in solid/solid two-phase systems [54]. This equation is derived from Darken's method [64], which can be found in [65], and is based on the intrinsic fluxes of the atoms and a drift velocity for mass conservation [66–69]. Its general formulation is:
$$
\left\{
\begin{array}{l}
\displaystyle \frac{\partial c_k}{\partial t} + V \cdot \nabla c_k = -\nabla \cdot (\mathbf{j}_k^d) - c_k \nabla \cdot V \quad (a) \\
\\
\qquad \text{with} \\
\\
\displaystyle \mathbf{j}_k^d = -c_k \frac{D_k}{RT} \nabla \mu_k \quad (b) \\
\\
\displaystyle \mu_k(T, p, x_k) = \mu_k^*(T, p^*, x_k^*) + RT \ln\left(\frac{x_k}{x_k^*}\right) + \Omega_k(p - p^*) \quad (c) \\
\\
\displaystyle V = \mathbf{V}^D + \mathbf{V}^\sigma \quad (d) \\
\\
\displaystyle \mathbf{V}^D = -\sum_{k=A,B} \Omega_k \mathbf{j}_k^d \quad (e)
\end{array}
\right.
\quad (4)
$$

where: i) $c_k$ is the concentration of the k-species; ii) $\mathbf{j}_k^d$ is the intrinsic diffusion flux of the k-species defined using the intrinsic diffusion coefficient $D_k$ and the extended chemical potential $\mu_k$; iii) $\mu_k$ is written using the pure chemical term $RT \ln(x_k/x_k^*)$ ($x_k$ and $x_k^*$ being the molar fraction and equilibrium molar fraction of the k-species, respectively) and the mechanical term $\Omega_k(p-p^*)$ ($p$ and $p^*$ being the pressure and reference pressure, respectively) and $\Omega_k$ the molar volume of the k-species). Note that $\mu_k^*(T, p^*, x_k^*)$ is the reference chemical potential which is equal for both phases at pressure $p^*$ and temperature $T$. Interface equilibrium implies that $\mu_{k\Sigma}^\alpha(T, p) = \mu_{k\Sigma}^\beta(T, p)$, where $\mu_{k\Sigma}^\alpha(T, p)$ and $\mu_{k\Sigma}^\beta(T, p)$ are the mechano-chemical potentials of the k-species at the interface on the $\alpha$-phase (solution) and the $\beta$-phase (solid) sides, respectively. Assumptions are made about the molar volumes of the species, which are considered to be identical for all species in the two phases ($\Omega_A^\alpha = \Omega_B^\alpha = \Omega_B^\beta = \Omega$) and independent of pressure. This assumption is valid for condensed phases, which is the case in the framework of this article, which deals with liquids and solids and implies that the concentrations of species are not affected by variations in pressure. In that case, the k-species molar fractions $x_k$ and $x_k^*$ used to define the chemical potential $\mu_k$ can be replaced by the species concentrations $c_k$ and $c_k^*$. Then, by developing the mechano-chemical potential $\mu_k$, the $\beta$-phase being a pure phase, we found at the interface:

$$
RT \ln \left( \frac{c_{k\Sigma}^{\alpha}}{c_{k}^{\alpha*}} \right) + \Omega(p_{\Sigma}^{\alpha} - p^{*}) = \Omega(p_{\Sigma}^{\beta} - p^{*}), \quad (5)
$$

which leads to Eq. (2) defined in Section 4.1 when $p_{\Sigma}^{\alpha} = p^{*}$.

iv) **V** is the advection velocity equals for every species and corresponding to the addition of the Darken velocity **V**^D and mechanical velocity **V**^σ. **V**^σ is either the velocity of elastic deformation of the solid or the viscous flow velocity of the liquid. The Darken velocity **V**^D is calculated with the intrinsic flux of each species and their molar volume Ω_k (Eq. (4e)). The velocity **V**^σ is the result of the mechanical governing equation solving. Refer to the literature for a full description of the unified formulation of an elastic solid and a compressive fluid [61].

The interface displacement is determined using the advection equation of the phase function $\phi$, it is written as:

$$
\frac{\partial \phi}{\partial t} = (\mathbf{V}^D + \mathbf{V}^\sigma) \cdot \nabla \phi = 0. \quad (6)
$$

Finally, for solving this diffusion transport equation, Neumann conditions are considered at system boundaries ensuring zero flux at boundaries and thus mass conservation.

### 4.2.3. Numerical methodology and discretization

The mechanical and diffusion equations are solved using an Eulerian method in the framework of fixed grid, in which the algebraic Volume-Of-Fluid (VOF)-based interface advection method [70] is used to differentiate the solid phase from liquid phase and to deal with interface displacement. The one-field equations are discretized in time and space by explicit/implicit methods, respectively, on a staggered grid [62,63]

using the parallel framework of the Notus Computational Fluids Dynamics code [71–74]. A second order centered scheme is used for the discretization of the diffusion term in Eq. (4a). A Lax-Wendroff scheme with a super bee flux limiter [70] is used to solve the advection term in Eqs. (4a) and (6).

## 4.3. Results

### 4.3.1. Validation of the mechano-chemical coupling at interfaces

In the first step, the dissolution/transport/precipitation mechanism is simulated by imposing pressure jumps on two solid/liquid interfaces, that remain constant during the simulation of the process. The initial molar fraction of B-species in the solution is that of saturated solution at pressure $p^*$ and temperature $T$, $x_B^{\alpha 0} = x_B^{\alpha}(p^*) = 0.5$. The molar fraction in solid 1 and solid 2 are $x_1^\beta = x_2^\beta = 1$. The initial molar fraction profile, when all the phases are under pressure $p^*$, is represented in Fig. 8b by the dotted black line. When the pressure jumps ($p_{\Sigma1}^\beta - p^*, i = 1, 2$) are imposed at the interfaces $i = 1, 2$, the equilibrium molar fractions of the solute are shifted towards higher values, with $x_{\Sigma1}^\alpha > x_{\Sigma2}^\alpha$ in accordance to $p_{\Sigma1}^\beta > p_{\Sigma2}^\beta$. The map of the molar fraction of species B in the solution is shown at different characteristic times during the process (Fig. 8a). It is worth noting that these molar fractions at the interfaces induced by mechano-chemical equilibrium remain constant during the process, as shown by the concentration profiles (Fig. 8b). Their values correspond perfectly to the analytical values calculated from Eq. (2) and shown in Table 3. Note that when total species concentration $c$ is constant for any points of the system and for any pressure (valid for condensed matter), the concentrations in Eq. (2) can be replaced by molar fractions ($c_{\Sigma1}^\alpha = c_{\Sigma2}^\alpha$ and $c_{\Sigma}^\alpha = c_{\Sigma}^\alpha$).

The simulated dissolution-precipitation process highlights 3 stages: i) the first one for which the solution is undersaturated relative to the pressure induced equilibrium mole fractions $x_{\Sigma1}^\alpha$ and $x_{\Sigma2}^\alpha$ of the two interfaces (red curve in Fig. 8b). The solid phases therefore dissolve at both interfaces and the volume of the liquid phase increases; ii) the second in which the solution is saturated near the interface of solid 2 due to a zero molar fraction gradient. Interface 2 is thus immobile under these conditions (green curve in Fig. 8b); iii) from this point onwards, due to the transport of species B towards interface 2, the solution near this interface 2 becomes supersaturated in relation to the equilibrium molar fraction $x_{\Sigma2}^\alpha$. A precipitation process takes place at this interface 2, while dissolution continues to occur at interface 1 (purple curve in Fig. 8b). In this last stage the two interfaces move in the same direction.

### 4.3.2. Simulation of mechano-chemical coupling for a porous medium subjected to uniaxial pressing

In this section, we aim to simulate more realistic dissolution-precipitation processes by no longer considering fixed pressure jumps at the liquid/solid interfaces, but pressure jumps derived from stress field calculations in a porous elastic solid subjected to uniaxial pressure. Fig. 9a shows the pseudo 3D configuration of the porous solid characterized by pores aligned in a direction perpendicular to the axis of pressure application. These pores are filled with liquid and the pressure in the liquid is kept constant at the reference pressure $p^*$.

This 2D-configuration can be viewed as open porosity segments oriented perpendicularly to the uniaxial applied pressure. The solid phase is considered to be a continuous medium with no intergranular liquid film. Only the pore surface is involved in the pressure solution mechanism. The relative density of 77 % is close to the experimental one (i.e. 73 %). For this relative density, the pore surface of silica is expected to be totally wetted with liquid layer thick enough for the silica solute to diffuse through it from one to another point located on the solid/liquid surface (see thermodynamic results in Section 3.2). So, the condition of a fully saturated porosity also allows the solute to diffuse easily to reach any point on the surface. The constant liquid pressure $p^*$, is
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** 1D-simulation of the dissolution/transport/precipitation process. The middle phase of the system is a solution of B-species at pressure $p^*$, which is between pure solids (orange color) of B-species. The blue scale corresponds to the molar fraction of the solute in the liquid solution. At initial time, the initial molar fraction of the solute is that of the saturated solution at $p^*$, $x_B^{t0} = x_B^{t*}(p^*) = 0.5$. The phases located at the two ends correspond to pure solid phases $x_1^β = x_2^β = 1$, whose respective pressures have been fixed and kept constant during the process evolution. The values of the pressure differences considered between the liquid and the solid are $p_{\Sigma1}^β - p^* = 70$ MPa (solid on the right), and $p_{\Sigma2}^β - p^* = 20$ MPa (solid on the left). a) Characteristic configurations during the evolution of the process at temperature $T = 300$ °C showing the two interface displacements and the molar fraction map; b) corresponding molar fraction profiles at same non-dimensional evolution times.</figcaption></figure>

<table><caption>**Table 3**

Mechano-chemical equilibrium of the interface. The molar fraction displacements of the solute at interface induced by pressure jumps are calculated from Eq. (2) by considering $T = 300$ °C and molar volume $\Omega = 2.0 \times 10^{-5} \text{m}^3/\text{mol}$.</caption><thead><tr><th></th><th>Interface 1</th><th>Interface 2</th></tr></thead><tbody><tr><td>Interface pressure jump (MPa)</td><td>70</td><td>20</td></tr><tr><td>Displaced molar fraction</td><td>0.67</td><td>0.54</td></tr></tbody></table>

representative of the saturation vapor pressure $p_S^T$ at temperature $T$, which also remains constant, whatever the volume of the pores for the two-phase pore space (refer to Section 3.2 concerning the pore space thermodynamics).

The elastic constants of the solid and the compressibility of the fluid are provided in the legend to Fig. 9, as well as the dynamic constants, diffusion coefficients for the solid and liquid and viscosity of the liquid. Note that the value of the isothermal compressibility coefficient $\chi_T^T$ of the liquid is very low, characterizing thus a fluid close to incompressibility. The displacements of the pore surface, induced by mass transfers, lead to liquid flow. The non-dimensional ratio $\frac{\rho_D^{D'}}{\rho^0} = 5 \times 10^{-4}$ considered, means that the interface displacement which is limited by the diffusion of species in the liquid, is also the kinetically limiting process with respect to the liquid flow. These conditions mean that during the process, the liquid is in mechanical equilibrium, i.e. at constant pressure. As

the volume of the pores remains constant throughout the process, the pressure does not change either.

Fig. 9b shows the pressure field $(p^β - p^*)$ calculated in the solid phase and resulting from the uniaxial compression of the porous solid. For the feasibility of the simulation, the pressure range in the solid is underestimated compared to that expected for uniaxial pressing at 350 MPa, and the value of the molar fraction of the silica solute in saturated solution is overestimated. Nevertheless, these configurations would give the same trend for the microstructural evolution, only the kinetics would be affected. As expected, this pressure map shows an anisotropic stress field in the solid which induces pressure jump gradients at the solid/liquid interfaces. These pressure jumps at the interfaces are higher on the horizontal lines through the pores than on the vertical lines. In the light of the results reported in Section 4.3.1 and the characteristics of the pressure jumps along the interfaces, the mechano-chemical equilibrium solute molar fraction should be more strongly displaced on the interface points located on the horizontal lines passing through the pores than on the points located on the vertical lines (see Fig. 9c). Thus dissolution/transport/precipitation mechanism should be activated (Fig. 9c).

As expected, on the basis of the elastic stress map of the solid in the initial state, the simulation of mass transport and transfer at solid/liquid interfaces shows dissolution of the solid where the pressure jumps at the interface are highest, and precipitation where the pressure jumps are lowest (Fig. 10a). Stress anisotropy leads to an anisotropic morphological evolution of the pores, which can be clearly seen by superimposing
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** a) Initial configuration of the porous elastic solid with pore space filled with liquid water; b) Initial pressure map ($p-p^*$) of the elastic solid uniaxially pressed. Elastic constants considered for the solid are: Young's modulus $E_V^\beta = 7 \times 10^{10}$ Pa; Poisson coefficient $\nu = 0.17$, 1st Lamé's coefficient $\mu_E^\beta = 1.5 \times 10^{-10}$ Pa, isothermal compressibility coefficient $\chi_T^\beta = 2.8 \times 10^{-11}$ Pa$^{-1}$ and for the liquid $\chi_i^\beta = 3 \times 10^{-9}$ Pa$^{-1}$; c) expected dissolution/transport/precipitation process in reference to mechano-chemical coupling at the interface and related Eq.(2). The value of the solute molar fraction of the saturated solution at $p^*$ is $x_B^{\beta*}(p^*) = 0.5$, and the initial molar fraction is $x_B^0 = 0.51$; the dynamic constants normalized to the system dimension $L$ are for the solid the diffusion coefficient $D^a/L^2 = 5 \times 10^{-4}$ s$^{-1}$ and for the liquid $\eta^a/(\rho L^2) = 1.0$ s$^{-1}$.</figcaption></figure>

the pore surface plot in the initial state (green line). This change in pore shape occurs with a variation in pore volume of less than 2% (Fig. 10b). The equilibrium molar fractions $x_1^\alpha$ and $x_2^\alpha$ for the largest ($p_{\Sigma1}^\beta - p^*$) and smallest ($p_{\Sigma2}^\beta - p^*$) pressure jumps are calculated from Eq. (2). It is found that the initial mole fraction is not centred with respect to these values, $x_{B\Sigma2}^\alpha = 0.501 < x_B^0 = 0.51 < x_{B\Sigma2}^\beta = 0.54$. This means that undersaturation is greater than oversaturation. In other words, and in these conditions, dissolution of the solid will predominate over precipitation. This explains the pore volume increase at the beginning of the process as it was the case for the validation test (Section 4.3.1). After this transition period, pore volume remains almost constant. The anisotropic evolution of pore morphology at constant volume implies that they can be characterised by two lengths, one larger than the initial diameter and the other smaller.

## 5. Discussion on the relevance of the 2D simulation to the evolution of the pore size distribution and its correlation with experimental results

In this section, we discuss all the experimental, modelling and simulation results obtained in order to support our arguments that the 2D simulation provide reliable insight into the mechanisms of silica hydrothermal sintering.

### 5.1. On the question of liquid water distribution

As pointed out in the introduction, the solvent addition is the key point towards lowering sintering temperatures by activating dissolution-precipitation processes, either for CSP or HS sintering. As HS evolves in a closed system, it was possible to quantify the degree of liquid water saturation in the porosity using the thermodynamic properties of water (see Section 3). The modelling diagram in Fig. 6, clearly shows that for 7

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Evolutions of the pore morphology and pore volume of a porous medium uniaxially pressed. a) pressure map ($p-p^*$) at dimensionless time $t \cdot D/L^2 = 0.015$ of the evolving process; green line corresponds to the initial solid/solution interface and black line to the solid/solution interface at time $t \cdot D/L^2 = 0.015$. b) evolution of the total pore volume $V_P$ normalised to the initial pore volume $V_P^0$.</figcaption></figure>
– 10 wt. % water present initially in the powder (experimentally evaluated), the liquid water is not totally evaporated at T = 300 °C. More precisely, it has been estimated that around 50 % of the water introduced remains in a liquid state at this temperature, with a limited relative density of 75 %. Since the chemical interactions between silica and water favour silica wetting, it's reasonable to assume that the remaining liquid water covers the entire pore surface with sufficient layer thickness for the solute to diffuse from any point on the pore surface to any other on that surface.

## 5.2. On the question of local solute transport

In this likely hypothesis of an entirely wetted pore surface deduced from thermodynamic modelling, the junction between solid, liquid and vapor phases is absent, preventing thus the capillary flow along the connected porosity. In addition, as the total porosity is not fully saturated with liquid water, any reduction in pore volume does not induce long-distance flow through the porosity either. These two last points mean that long-distance solute transport driven by liquid flow is limited. Only solute transport by diffusion driven by concentration gradients in the solution should predominate, which implies solute transport over short distances. Such conditions strongly suggest that during the densification plateau observed at 300 °C (see Section 2.2 Experimental protocol - Fig. 1), only local transport is involved.

## 5.3. On the question of porous microstructure configuration

The relative density around 75 % corresponds to the intermediate stage of sintering, which is mainly characterized by interconnected pores [75]. In Coble's geometric model, this connected porosity is represented by channels located on the edges of polyhedra, bridging larger pores on the vertices. Microtomography images obtained during solid-state sintering of spherical soda-lime glass particles with a size range of 100–200 µm [76,77], exhibited such channels at a density of around 72 %. As the material is amorphous, these channels were observed to evolve into cylindrical shape before breaking down for close pore formation at a larger relative density above 85 %. In the simple Coble's representation also considering cylindrical connected channels [75], some of them are oriented in the same direction. The 2D configuration of aligned cylindrical pores tends to represent this category of channels present in the simple representation of connected porosity.

## 5.4. On the link between pore morphology and pore size distribution

The results obtained from mercury intrusion porosimetry highlighted a change in pore morphology evolution without densifying mechanism (Section 2.2). Pore size distribution is calculated by applying the Washburn equation, assuming that pores are cylinders of different diameters. When the pores are no more cylindrical, they can be characterized by two diameters, a small one defined as the minimum distance between two distinct surface points with identical tangents, and a large one defined similarly as the maximun distance. A reduction in the small diameter tends to increase the mercury inlet pressure. The 2D simulation of the morphology evolution of cylindric pores (Fig. 10), which mimics some channels in the connected porosity subjected to uniaxial compression, clearly shows such small diameter reduction compatible with a pressure inlet increase. When mercury enters the interconnected open porosity, it necessarily passes through these specific channels to fill the entire pore volume like in a “bottle-neck” pore configuration. This local simulation of the evolution of pore morphology is consistent with the experimental result of an increase in the proportion of mesopores at constant relative density.

## 5.5. General comments on the HS of silica

Pure water is not an efficient solvent for dissolving silica, which does

not facilitate the diffusion and the diffusion/sliding mechanisms at grain boundaries inherent to the densification process. Even if large stress gradients are present at the grain boundaries, diffusion coefficients can become very low due to the limited penetration of water between particles. Transport at the grain boundaries then becomes very limited, which contributes to the mechanical stabilisation of the porous structure and the annihilation of densification. Nevertheless, when a layer of water is present on the overall pore surface, the dissolution/precipitation mechanism via diffusion of species in the liquid phase remains possible, even though the stress gradients at the pore surface are much lower than at the contacts between the grains. This low stress gradient is balanced by a diffusion coefficient that can allow matter redistribution in porosity at constant pore volume.

# 6. Conclusion

The case study performed on hydrothermal sintering of silica powder with water as a solvent showed results consistent with those reported in the literature, in particular in terms of the relative density obtained that do not exceed 75 %. Here, considering a grain size of 300 nm diameter, the addition of 7–10 wt. % water and an applied uniaxial pressure of 350 MPa, the relative density of the silica powder increases from 52 % to 73 ± 2 % during the heating ramp until the dwell temperature of 300 °C is reached. The relative density then remains constant during the dwell temperature. Thermodynamic modelling of the pore space showed that, for the conditions of hydration of the green compacts, the degree of liquid water saturation in porosity remained high (28 – 50 %) at T = 300 °C and d = 75% due to the closed-system sintering conditions that limit evaporation. The observed stagnation in the densification is therefore not linked to an absence of liquid water in the compact, but rather to a very low hydration of the inter-granular contact limiting the pressure solution creep mechanism in this area.

However, an evolution of the pore size distribution was highlighted through mercury intrusion porosimetry during this plateau of relative density. Based on these macroscopic observations, we therefore assumed a change in the morphology of the pores at a constant pore volume. The scale of the pore was the focus of our study, which involved to develop specifically an original methodology for direct numerical modelling of the dissolution/precipitation mechanisms based on a mechano-chemical coupling of the transfer at the solid/liquid interfaces. Note that, to the best of our knowledge, it the first time that such direct simulations of mass transfer at the interfaces are applied to hydrothermal sintering. If the 2D numerical modelling is inherently limited in its ability to describe 3D microstructure evolutions, we demonstrate that the numerical simulation of the pressure-induced matter redistribution along the pore surface of a 2D porous medium fully saturated with liquid and uniaxially pressed was nevertheless relevant to interpret the experimental observations. The simulated morphology changes of the local pore sections were resulting in narrowing the pore size in some directions and widening in others ones. Such morphology changes were shown to be consistent with an artificial increase of the mesoporous proportion while no densification took place as it was observed experimentally.

To sum up, after a partial densification, the HS of silica with 7–10 wt. % pure water is most likely dominated by local dissolution/precipitation mechanism at the solid/liquid interface of the pores. The identification of such non-densifying mechanism, through a simulation/experimental coupling strategy, represents a first in the field of HS.

## CRediT authorship contribution statement

S. Bordère: Writing – review & editing, Writing – original draft, Supervision, Software, Methodology, Conceptualization. L. Villatte: Writing – review & editing, Visualization, Validation, Investigation. M.-A. Dourges: Writing – review & editing, Visualization, Validation, Investigation. G. Goglio: Writing – review & editing, Project administration, Funding acquisition. C. Elissalde: Writing – review & editing,
Writing – original draft, Investigation, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work has been made possible by the financial support of the French National Research Agency (ANR) in the framework of the project HYDILIC (ANR-18-CE08-0006). Numerical simulations were carried out using the facilities of the MCIA CURTA mesocentre.

Appendix A. Nomenclature

A.1. Sintering parameters

$T_0$: room temperature

$T$: sintering dwell temperature

$d_0$: initial powder relative density

$d$: powder relative density

A.2. Mass

$m_S$: mass of the solid into cell (invariant)

$m_L^{T_0}$: mass of the liquid solvent into cell at room temperature $T_0$ (invariant)

$m$: total mass of the solvent (liquid and vapor) distributed into porosity and reservoir volumes at (T, $p_A$)-conditions; $m = m_L^{T_0}$ (invariant)

A.3. Characteristic volumes

$V_S$: volume of the solid phase (invariant)

$V_R$: volume of reservoirs (invariant)

$V_L^{T_0}$: volume of the liquid solvent at room temperature $T_0$

$V_C(d)$: volume of the sampling cell at relative density $d$

$V_P^T(d)$: volume of the porosity at relative density $d$

$V_F(d)$: volume of the fluid phases (liquid and vapor) including reservoir volume at relative density $d$

$V_L^T(d)$: volume of the liquid solvent at room temperature $T$ and relative density $d$

A.4. Densities

$\rho_L^{T_0}$: density of the liquid solvent at room temperature $T_0$

$\rho_V^{T_0}$: density of the vapor solvent at room temperature $T_0$

$\rho_S^{T_0}$: density of the solid at room temperature $T_0$

$\rho_M^T(d)$: apparent density of the solvent (liquid and vapor) relative to porosity and reservoir volumes at temperature $T$ and powder relative density $d$;

$\rho_M^T(d) = m/V_F^T(d)$

$\rho_L^T$: density of the solvent in liquid phase at temperature $T$ for thermodynamic saturation points

$\rho_V^T$: density of the solvent in vapor phase at temperature $T$ for thermodynamic saturation points

$\rho_S^T$: density of the solid at temperature $T$; $\rho_S^T \approx \rho_S^{T_0}$ (invariant)

A.5. Mass fraction

$f_{mL}^{T_0}$: Initial mass fraction of the liquid solvent in cell at room temperature $T_0$

A.6. Saturation of the porosity volume with liquid solvent

$f_L^{T_0}(d_0)$: Initial volume fraction of the liquid solvent in porosity at room temperature $T_0$ and relative density $d_0$

$f_L^T(d)$: Volume fraction of the liquid solvent in porosity at temperature $T$ and relative density $d$

A.7. Saturation of the fluid volume with liquid solvent

$x_L^{T_0}(d_0)$: Initial volume fraction of the liquid solvent in the overall fluid volume $V_F(d_0)$ at room temperature $T_0$ and relative density $d_0$

$x_L^T(d)$: Volume fraction of the liquid solvent in the overall fluid volume $V_F(d)$ at temperature $T$ and relative density $d$

Appendix B. Calculation of the liquid solvent saturation of porosity at sintering temperature T and density d

The variables and constants of the problem are defined in Appendix A.

First, note that the saturation of volume porosity with solvent at temperature $T$ is independent of the applied pressure $p_A$. The calculation of this
saturation value only requires the knowledge of mass values $m_S$, $m_L^{T_0}$, powder relative densities values $d_0$, $d$, device volume values $V_R$, $V_C(d_0)$ and densities values $\rho_S^{T_0}$, $\rho_L^{T_0}$, $\rho_V^{T_0}$, and $\rho_S^T$, $\rho_L^T$, $\rho_V^T$. All these densities can be obtained from thermodynamic data of fluids at saturation points (equilibrium between the liquid and vapor phases of the solvent). Relationships for calculating the final saturation of volume porosity with liquid solvent are defined below.

## B.1. Calculation of the initial mass fraction of liquid solvent $f_{mL}^{T_0}$ at room temperature $T_0$

The mass fraction is defined from the masses of powder and added liquid solvent as:

$$f_{mL}^{T_0} = \frac{m_L^{T_0}}{m_L^{T_0} + m_S} \qquad (B.1)$$

## B.2 Calculation of the volume saturation of liquid solvent in pore space $f_L^{T_0}$ at room temperature $T_0$ and initial powder relative density $d_0$

At room temperature, for aqueous solutions, we have $\rho_V^{T_0} \ll \rho_L^{T_0}$. $\rho_V^{T_0}$ is thus neglected into the conversion relation of mass fraction of liquid solvent $f_{mL}^{T_0}$ to saturation volume $f_L^{T_0}$ of liquid solvent in pore space.

$$f_L^{T_0}(d_0) = \frac{f_{mL}^{T_0} d_0 \rho_S^{T_0}}{(1 - d_0)((1 - f_{mL}^{T_0})\rho_L^{T_0})} \qquad (B.2)$$

## B.3 Calculation of characteristic volumes related to the sintering device

From the relative density $d_0$, the volume of the sampling cell $V_C(d_0)$ and from the initial volume fraction of the liquid solvent in porosity $f_L^{T_0}(d_0)$, the volumes of the solid $V_S$ (invariant) and the liquid $V_L^{T_0}(d_0)$ are defined as,

$$V_S = d_0 V_C(d_0) = d V_C(d) \quad (\text{invariant}), \qquad (B.3)$$

$$V_L^{T_0}(d_0) = V_C(d_0)(1 - d_0)f_L^{T_0}(d_0). \qquad (B.4)$$

The porosity volume $V_P^T(d)$ at every powder relative density $d$ is,

$$V_P^T(d) = V_C^T(d)(1-d) = V_C(d_0) \frac{d_0}{d}(1-d), \qquad (B.5)$$

and the total fluid volume $V_F^T(d)$, by considering the volume of all reservoirs $V_R$ (invariant) is,

$$V_F^T(d) = V_C^T(d)(1-d) + V_R = V_C(d_0) \frac{d_0}{d}(1-d) + V_R. \qquad (B.6)$$

## B.4. Calculation of the apparent density of solvent $\rho_M^T(d)$ relative to porosity and reservoir volumes at temperature $T$ and powder density $d$

$$\rho_M^T(d) = \frac{m}{V_F^T(d)} = \frac{\rho_L^{T_0} V_L^{T_0}(d_0)}{V_F^T(d)}. \qquad (B.7)$$

## B.5. Calculation of the liquid solvent saturation relative to total fluid volume (porosity and reservoir volumes) and porosity volume

The liquid solvent saturations in both volumes are determined from the involved solvent thermodynamics. In the following, we detail the different calculation steps.

In Fig. B.1, is drawn using the REFPROP database [38] the pressure-density phase diagram of pure water as example. On this diagram, three isotherm lines are shown. From the knowledge of the apparent density $\rho_M^T(d)$ (Eq. (B.7)) at temperature $T$, we illustrate how we can determine whether the pore space is characterized by a single-phase or two-phase domain.

More precisely, if the apparent solvent density $\rho_M^T(d)$ is:

i) below the vapor phase density $\rho_V^T$, (for instance $\rho_M^T(d) = \rho_1^T$, see Fig. B.1) the pore space is monophasic with a vapor pressure $p_1^T$ depending on temperature $T$ and $\rho_M^T(d)$ (bivariant system). In that case, the values of the volume fractions of liquid solvent in the overall fluid volume and porosity volume are respectively, $x_L^T(d) = 0$ and $f_L^T(d) = 0$;

ii) similarly, above the liquid phase density (for instance $\rho_M^T(d) = \rho_2^T$, see Fig. B.1), the pore space is also monophasic with a liquid pressure $p_2^T$ depending on $T$ and $\rho_M^T(d)$ (bivariant system). In that case, $x_L^T(d) = 1$ and $f_L^T(d) = 1$;
iii) in-between the two densities of the vapor and liquid phases, respectively $\rho_V^T$ and $\rho_L^T$, (for instance $\rho_M^T(d) = \rho_D^T$, see Fig. B.1) the pore space is diphasic. Whatever the apparent solvent density in this domain, the pressure inside the porosity equals the saturation pressure $p_S^T$, that is the system only depends on the temperature $T$ (monovariant system). Within this two-phase domain, the liquid solvent saturation relative to the total fluid volume is defined as,

$$x_L^T(d) = \frac{(\rho_M^T(d) - \rho_V^T)}{(\rho_L^T - \rho_V^T)}, \qquad (B.8)$$

and that relative to the porosity volume as,

$$f_L^T(d) = x_L^T(d) \frac{V_F^T(d)}{V_P^T(d)}. \qquad (B.9)$$

During densification, as the fluid volume $V_F^T(d)$ decreases, the apparent solvent density $\rho_M^T(d)$ increases giving rise to mass transfer from vapor phase to liquid phase at constant pressure $p_S^T$. In other words, the liquid volume increases in a decreasing porosity volume, two evolutions that each contributes to increase the liquid saturation.

<figure><img src="image_13.png" /><figcaption>Fig. B.1. Pressure-density phase diagram of pure water [38]. Three isotherms are drawn on this diagram. The position on this diagram of the apparent solvent density $\rho_M^T(d)$ of the fluid occupying the pore and reservoir volumes can be used to determine the volume fraction properties of the pore space. For $T = 300$ °C, if $\rho_M^T(d) = \rho_1^T < \rho_V^T$, the pore space is a single vapor phase at pressure $p_1^T$; if $\rho_V^T < \rho_M^T(d) = \rho_D^T < \rho_L^T$ the pore space is biphasic at the saturation pressure $p_S^T$; if $\rho_M^T(d) = \rho_2^T > \rho_L^T$, the pore space is a liquid single-phase at pressure $p_2^T$.</figcaption></figure>

## Appendix C. Calculation of the thickness of the solvent layer

Let's consider a spherical morphology of particles of radius $r$ uniformly covered by a layer of liquid solvent of thickness $e$. Before densification and at any temperature, the ratio between the volume of liquid $V_L^{T_0}(d_0)$ and that of solid $V_S$ is approximated geometrically, when $e/r \ll 1$, by:

$$\frac{V_L^{T_0}(d_0)}{V_S} \approx \frac{3e}{r}. \qquad (C.1)$$

Replacing these volumes by their expression in the Eqs. (B.3) and (B.4), we obtained:

$$e \approx \frac{1 - d_0}{3d_0} f_L^T(d_0) r. \qquad (C.2)$$

Fig. C.1 shows the proportion of water layer thickness to particle radius that is numerically calculated for $f_L^{T_0}(d_0) = 0.40$ and $d_0 = 0.60$. This calculation takes account of the particle contact $x_r = 0.3$, which will slightly increase the liquid layer thickness compared with the analytical calculation (Eq. (C.2)). The value of 0.09 is obtained thus validating the analytical equation giving rise to $e/r = 0.088$.
<figure><img src="image_14.png" /><figcaption>Fig. C.1. Visualisation of the proportion of the liquid layer thickness to the particle radius numerically obtained for the two-particle configuration with a particle contact $x_r = 0.3$ considering for liquid saturation level of porosity the value $f_l^{T_0}(d_0) = 0.40$ at room temperature $T_0$ and relative density $d_0 = 0.60$.</figcaption></figure>

References

[1] L. Karacasulu, M. Karakava, U. Adem, V.M. Slgavo, M. Biesuz, C. Vakifahmetoglu, Fast-firing of potassium sodium niobate (KNN), Open Ceramics 17 (2024) 100541.

[2] J.L. Miranda, C. Ribeiro, A. Tkach, A.M.O.R. Senos, P.M. Vilarinho, Ultrafast high-temperature sintering of potassium-sodium niobate: processing and properties, J. Eur. Ceram. Soc. 44 (2024) 116701.

[3] O. Guillon, W. Rheinheimer, M. Bram, A perspective on emerging and future Sintering technologies of ceramic materials, Adv. Eng. Mater. 25 (2023) 2201870.

[4] M. Biesuz, S. Grasso, V.M. Sglavo, What's new in ceramics sintering? A short report on the latest trends and future prospects, Curr. Opin. Sol. Stat. Mater. Sci. 24 (2020) 100868.

[5] T. Ibn-Mohammed, C.A. Randall, K.B. Mustapha, J. Guo, J. Walker, S. Berbano, S. C. Koh, D. Wang, D.C. Sinclair, I.M. Reaney, Decarbonising Ceramic Manufacturing: a Techno-Economic Analysis of Energy Efficient Sintering Technologies in the Functional Materials Sector, J. Eur. Ceram. Soc. 39 (2019) 5213-5235.

[6] J. Guo, H. Guo, A.L. Baker, M.T. Lanagan, E.R. Kupp, G.L. Messing, C.A. Randall, Cold sintering: a paradigm shift for processing and integration of ceramics, Angew. Chem. Int. 55 (2016) 11457-11461.

[7] J. Guo, R. Floyd, S. Lowum, J.-P. Maria, T. Herisson de Beauvoir, J.-H. Seo, C. A. Randall, Cold Sintering: progress, Challenges, and Future Opportunities, Annu. Rev. Mater. Res. 49 (2019) 275-295.

[8] A. Ndayishimiye, S.H. Bang, C.J. Spiers, C.A. Randall, Reassessing cold sintering in the framework of pressure solution theory, J. Eur. Ceram. Soc. 43 (2023) 1-13.

[9] A. Ndayishimiye, M.Y. Sengul, T. Sada, S. Dursun, S.H. Bang, Z.A. Grady, K. Tsuji, S. Funahashi, A.C.T. van Duin, C.A. Randall, Roadmap for densification in cold sintering: chemical pathways, Open Ceramics 2 (2020) 100019.

[10] J. Guo, S.S. Berbano, H. Guo, A.L. Baker, M.T. Lanagan, C.A. Randall, Cold Sintering Process of Composites: bridging the Processing Temperature Gap of Ceramic and Polymer Materials, Adv. Funct. Mater. 26 (2016) 7115-7121.

[11] A. Jabr, H.N. Jones, A.P. Argüelles, S. Trolier-McKinstry, C. Randall, R. Bermejo, Scaling up the cold sintering process of ceramics, J. Eur. Ceram. Soc. 43 (2023) 5319-5329.

[12] E.T. Gecil, A.A. Raja, Densification of ceramics and ceramics-based composites using ultralow temperature sintering (cold sintering): a comprehensive review, Int J Appl Ceram Tech 22 (2025) 14996.

[13] M. Yoshimura, S. Somiya, Fabrication of Dense, Nonstabilized ZrO2 Ceramics by Hydrothermal Reaction Sintering, in: S. Somiya, Y. Moriyoshi (Eds.), Sintering Key Papers, Springer Netherlands, Dordrecht, 1990, pp. 663-665.

[14] G. Goglio, A. Largeteau, A. Ndayishimiye, M. Prakasam, Procédé et Dispositif de Densification de Matériaux Ou de Consolidation d'un Assemblage de Matériaux Par Frittage Hydrothermal Ou Solvothermal, FR17/54585. 2017.

[15] G. Goglio, A. Largeteau, A. Ndayishimiye, M. Prakasam, France Brevet FR 3066759. 2018, International Patent WO 2018215559 (2018-11-29).

[16] G. Goglio, A. Ndayishimiye, A. Largeteau, C. Elissalde, View point on hydrothermal sintering: main features, today's recent advances and tomorrow's promises, Scripta Mater. 158 (2019) 146-152.

[17] A. Ndayishimiye, A. Largeteau, M. Prakasam, S. Pechev, M.-A. Dourges, G. Goglio, Low temperature hydrothermal sintering process for the quasi-complete densification of nanometric α-quartz, Scr. Mater. 145 (2018) 118-121.

[18] G. Goglio, A. Ndayishimiye, C. Elissalde, C.A. Randall, Cold Sintering and Hydrothermal Sintering. Reference Module in Materials Science and Materials Engineering, Elsevier, 2020, pp. 311-326.

[19] A. Ndayishimiye, M.Y. Sengul, S.H. Bang, K. Tsuji, K. Takashima, T. Hérisson de Beauvoir, D. Denux, J.-M. Thibaud, A.C.T. van Duin, C. Elissalde, G. Goglio, C. A. Randall, Comparing hydrothermal sintering and cold sintering process:

mechanisms, microstructure, kinetics and chemistry, J. Eur. Ceram. Soc. 40 (2020) 1312-1324.

[20] M. Biesuz, G. Taveri, A.I. Duff, E. Olevsky, D. Zhu, C. Hu, S. Grasso, Adv. Appl. Ceram 110 (2020) 75-89.

[21] G. Taveri, S. Grasso, F. Gucci, J. Toušek, I. Dlouhy, Bio-Inspired Hydro-Pressure Consolidation of Silica, Adv. Funct. Mater. 28 (2018) 1805794.

[22] M.P.A. van den Ende, G. Marketos, A.R. Niemeijer, C.J. Spiers, Investigating Compaction by Intergranular Pressure Solution Using the Discrete Element Method, J. Geophys. Res. Solid 123 (2018) 107-124.

[23] J.-P. Gratier, D.K. Dysthe, F. Renard, Chapter 2 - The Role of Pressure Solution Creep in the Ductility of the Earth's Upper Crust, Adv. Geophys. 54 (2013) 47-179.

[24] S. de Meer, C.J. Spiers, Uniaxial compaction creep of wet gypsum aggregates, J. Geophys. Res. Solid Earth. 102 (1997) 875-891.

[25] A.R. Niemeijer, C.J. Spiers, B. Bos, Compaction creep of quartz sand at 400-600 °C: experimental evidence for dissolution-controlled pressure solution, Earth Planet. Sci. Lett. 195 (2002) 261-275.

[26] C.J. Spiers, S. De Meer, A.R. Niemeijer, X. Zhang, Kinetics of rock deformation by pressure solution and the role of the thin aqueous films, in: Physicochem. Water Geol. Biol. Syst., S. Nakasimha, C. J. Spiers, L. Mercury, P. A. Fenter, M. F. J Hochella, Tokyo, Japan, 2004: pp. 129-158.

[27] S.H. Bang, A. Ndayishimiye, C.A. Randall, Anisothermal densification kinetics of the cold sintering process below 150 °C, J. Mater. Chem. C 8 (2020) 5668-5672.

[28] N. Albar, C. Manière, T. Hérisson de Beauvoir, E. Sanchez, G. Chevallier, A. Weibel, C. Estournès, Acta Mater. 290 (2025) 120974.

[29] L. Villatte, M.-I. Rua-Taborda, A. Ndayishimiye, C.A. Randall, A. Largeteau, G. Goglio, C. Elissalde, S. Bordère, Mechanisms and energetics in the early stages of solvent-assisted low-temperature sintering of ZnO, Materialia 22 (2022) 101418.

[30] N. Yamasaki, K. Yanagisawa, M. Nishioka, S. Kanahara, A hydrothermal hot-pressing method: apparatus and application, J. Mater. Sci. Lett. 5 (1986) 355-356.

[31] K. Yanagisawa, M. Nishioka, K. Ioku, N. Yamasaki, Neck formation of spherical particles by hydrothermal hot pressing, J. Mater. Sci. Lett. 10 (1991) 7-8.

[32] K. Yanagisawa, N. Nishioka, K. Ioku, N. Yamasaki, Densification of silica gels by hydrothermal hot-pressing, J. Mater. Sci. Lett. 12 (1993) 1073-1075.

[33] A. Ndayishimiye, A. Largeteau, S. Mornet, M. Duttine, M.-A. Dourges, D. Denux, M. Verdier, M. Gouné, T. Hérisson de Beauvoir, C. Elissalde, G. Goglio, Hydrothermal Sintering for Densification of Silica. Evidence for the Role of Water, J. Eur. Ceram. Soc. 38 (2018) 1860-1870.

[34] L. Villatte, Frittage Hydrothermal Un Procédé De Frittage Basse température, Doctoral dissertation, Université de Bordeaux, Bordeaux, France, 2022.

[35] A. Ndayishimiye, Sur Un Nouveau Procédé De Frittage De Céramiques à Basse Température : Le Frittage hydrothermal. Développement et Approche mécanistique, Doctoral dissertation, Université de Bordeaux, Bordeaux, France, 2017.

[36] S.-J.L. Kang, Sintering: Densification, Grain Growth and Microstructure, Elsevier Butterworth-Heinemann, Amsterdam Boston London, 2005.

[37] D. Bernache-Assollant, Chimie-physique Du Frittage, Hermès, Paris (France), Paris, 1993.

[38] E.W. Lemmon, I.H. Bell, M.L. Huber, M.O. McLinden, NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP, Version 10.0, National Institute of Standards and Technology, Standard Reference Data Program, Gaithersburg, 2018.

[39] D.B. Asay, S.H. Kim, Evolution of the Adsorbed Water Layer Structure on Silicon Oxide at Room Temperature, J. Phys. Chem. B 109 (2005) 16760-16763.

[40] M.S. Paterson, Nonhydrostatic thermodynamics and its geologic applications, Rev. Geophys. 11 (1973) 355-389.

[41] R.B. De Boer, On the thermodynamics of pressure solution—Interaction between chemical and mechanical forces, Geochim. Cosmochim. Acta. 41 (1977) 249-256.
[42] B. Malvoisin, L.P. Baumgartner, Mineral Dissolution and Precipitation Under Stress: model Formulation and Application to Metamorphic Reactions, Geochem. Geophys. Geosystems. 22 (2021) e2021GC009633.

[43] Z. Liu, C. Zhou, B. Li, L. Zhang, Y. Liang, Effects of grain dissolution-diffusion sliding and hydro-mechanical interaction on the creep deformation of soft rocks, Acta Geotech. 15 (2020) 1219-1229.

[44] S. Fleckenstein, D.A. Bothe, Volume-of-fluid-based numerical method for multi-component mass transfer with local volume changes, J. Comput. Phys. 301 (2015) 35-58.

[45] G. Gennari, R. Jefferson-Loveday, S.J. Pickering, A phase-change model for diffusion-driven mass transfer problems in incompressible two-phase flows, Chem. Eng. Sci. 259 (2022) 117791.

[46] A. Schulz, C. Wecker, V. Inguva, A.S. Lopatin, E.Y. Kenig, A PLIC-based method for species mass transfer at free fluid interfaces, Chem. Eng. Sci. 251 (2022) 117357.

[47] J. Maes, C. Soulaine, A unified single-field Volume-of-Fluid-based formulation for multi-component interfacial transfer with local volume changes, J. Comput. Phys. 402 (2020) 109024.

[48] C.P. Zanutto, E.E. Paladino, F. Evrard, B. van Wachem, F. Denner, Modeling of interfacial mass transfer based on a single-field formulation and an algebraic VOF method considering non-isothermal systems and large volume changes, Chem. Eng. Sci. 247 (2022) 116855.

[49] K. Ammar, B. Appolaire, G. Cailletaud, F. Feyel, S. Forest, Finite element formulation of a phase field model based on the concept of generalized stresses, Comput. Mater. Sci. 45 (2009) 800-805.

[50] R. Duddu, D.L. Chopp, P. Voorhees, B. Moran, Diffusional evolution of precipitates in elastic media using the extended finite element and the level set methods, J. Comput. Phys. 230 (2011) 1249-1264.

[51] M. Cottura, B. Appolaire, A. Finel, Y.Le Bouar, Phase field study of acicular growth: role of elasticity in Widmanstätten structure, Acta Mater. 72 (2014) 200-210.

[52] A. Durga, P. Wollants, N.A. Moelans, Quantitative phase-field model for two-phase elastically inhomogeneous systems, Comput Mater. Sci. 99 (2015) 81-95.

[53] R.D. Kamachali, C. Schwarze, M. Lin, M. Diehl, P. Shanthraj, U. Prahl, I. Steinbach, D. Raabe, Numerical Benchmark of Phase-Field Simulations with Elastic Strains: precipitation in the Presence of Chemo-Mechanical Coupling, Comput. Mater. Sci. 155 (2018) 541-553.

[54] S. Bordère, H.-P. Van Landghem, A. Redjaïmia, M. Gouné, Morphological and Dimensional Evolution of Nanosized Amorphous Silicon Nitride in α-Fe: diffusional and Elastic Effects, Metals (Basel) 14 (2024) 442.

[55] Z. Xu, P. Meakin, Phase-field modeling of solute precipitation and dissolution, J. Chem. Phys. 129 (2008) 014705.

[56] J.P. Pereira Nunes, M.J. Blunt, B. Bijeljic, Pore-scale simulation of carbonate dissolution in micro-CT images, J. Geophys. Res. 121 (2016) 558-576.

[57] Z. Liu, Q. Xu, J. Yang, K.H. Luo, L. Shi, Microcontinuum approach to multiscale modeling of multiphase reactive flow during mineral dissolution, Phys. Rev. Fluids 9 (2024) 043801.

[58] S.G. Kim, Large-scale three-dimensional simulation of Ostwald ripening, Acta Mater. 55 (2007) 6513-6525.

[59] H. Ravasha, L. Vanherpeb, J. Vleugelsa, N. Moelansa, Three-dimensional phase-field study of grain coarsening and grainshape accommodation in the final stage of liquid-phase sintering, J. Eur. Ceram. Soc. 37 (2017) 2265-2275.

[60] A. Sac-Morane, M. Veveakis, H. Rattez, A Phase-Field Discrete Element Method to study chemo-mechanical coupling in granular materials, Comput. Methods Appl. Mech. Engrg. 424 (2024) 116900.

[61] S. Bordère, J.-P. Caltagirone, A unifying model for fluid flow and elastic solid deformation: a novel approach for fluid-structure interaction, J. Fluids. Struct. 51 (2014) 344-353.

[62] S. Bordère, J.-P. Caltagirone, Multi-physics and multi-time scale approach for modeling fluid-solid interaction and heat transfer, Comput. Struct. 164 (2016) 38-52.

[63] S. Bordère, S. Glockner, Numerical modelling of diffusion-controlled phase transformation using the Darken method: application to the dissolution/precipitation processes in materials, Comput. Mater. Sci. 186 (2021) 109944.

[64] L. Darken, Diffusion, Mobility and Their Interrelation through Free Energy in Binary Metallic Systems, Trans, AIME 175 (1948) 184-201.

[65] S.A. Sridhar, Commentary on "Diffusion, Mobility and Their Interrelation through Free Energy in Binary Metallic Systems, L.S. Darken: Trans. AIME 175 (1948) 543-562, 184ffMetall. and Mat Trans A 412010.

[66] J. Herrnring, B. Sundman, B. Klusemann, Diffusion-driven microstructure evolution in OpenCalphad, Comput. Mater. Sci. 175 (2020) 109236.

[67] J. Ågren, Diffusion in phases with several components and sublattices, J. Phys. Chem. Solids 43 (1982) 421-430.

[68] J. Andersson, J. Ågren, Models for numerical treatment of multicomponent diffusion in simple phases, J. Appl. Phys. 72 (1992) 1350-1355.

[69] M. Danielewski, B. Wierzba, Thermodynamically consistent bi-velocity mass transport phenomenology, Acta Mater. 58 (2010) 6717-6727.

[70] C. Hirsch, Numerical Computation of Internal and External Flows volume 2. Computational Methods for Inviscid and Viscous Flows, 1st ed., John Wiley & Sons, Chichester, England, 1990.

[71] Notus CFD Code, available on line: https://notus-cfd.org. (accessed on 2 January 2016).

[72] J. Picot, S. Glockner, Reduction of the discretization stencil of direct forcing immersed boundary methods on rectangular cells: the ghost node shifting method, J. Comput. Phys. 364 (2018) 8-48.

[73] A.M.D. Jost, S. Glockner, A. Erriguible, Direct numerical simulations of fluids mixing above mixture critical point, J. Supercrit. Fluids 165 (2020) 104939.

[74] S. Glockner, A.M.D. Jost, A. Erriguible, Advanced petascale simulations of the scaling up of mixing limited flow processes for materials synthesis, Chem. Eng. J. 431 (2022) 133647.

[75] R.L. Coble, Sintering of crystalline solids. I. Intermediate and final state diffusion models, J. Appl. Phys. 32 (1961) 789-792.

[76] D. Bernard, D. Gendron, J.-M. Heintz, S. Bordère, J. Etourneau, First direct 3D visualisation of microstructural evolutions during sintering through X-ray computed microtomography, Acta Mater. 53 (2005) 121-128.

[77] D. Gendron, Etude Numérique Et Expérimentale Du Frittage à L'échelle Du grain, Doctoral dissertation, Université de Bordeaux, Bordeaux, France, 2001.