J Antoon Van Hooft

Stéphane Popinet

Bas J H Van De Wiel

Antoon Van Hooft
email: j.a.vanhooft@tudelft.nl

Adaptive Cartesian Meshes for Atmospheric Single-Column Models, a study using Basilisk 18-02-16

It is well known that the representation of certain atmospheric conditions in climate and weather models can still suffer from the limited grid resolution that is facilitated by modern-day computer systems. Herein we study a simple onedimensional analogy to those models by using a Single-Column Model description of the atmosphere. The model employs an adaptive Cartesian mesh that applies a high-resolution mesh only when and where it is required. The so-called adaptive-grid model is described and we report on our findings obtained for tests to evaluate the representation of the atmospheric boundary layer, based on the first two GABLS intercomparison cases. The analysis shows that the adaptive-grid algorithm is indeed able to dynamically coarsen and refine the numerical grid whilst maintaining an accurate solution. This is an interesting result as in reality, transitional dynamics (e.g. due to the diurnal cycle or due to changing synoptic conditions) are the rule rather than the exception.

Introduction

Single-Column Models (SCMs) are often used as the building blocks for Global (or General) Circulation Models (GCMs).

As such, many of the lessons learned from SCM development can be inherited by GCMs and hence the evaluations of SCMs receive considerable attention by the geoscientific model development community (see e.g. [START_REF] Neggers | Continuous single-column model evaluation at a permanent meteorological supersite[END_REF][START_REF] Bosveld | The third GABLS intercomparison case for evaluation studies of boundary-layer models. Part B: results and process understanding[END_REF][START_REF] Baas | From Near-Neutral to Strongly Stratified: Adequately Modelling the Clear-Sky Nocturnal Boundary Layer at Cabauw, Boundary-Layer Meteorology[END_REF]. In this work, we present a SCM that employs an adaptive Cartesian mesh that can drastically reduce the computational costs of such models, especially when pushing the model's resolution. The philosophy is inspired by recently obtained results on the evolution of atmospheric turbulence in a daytime boundary layer using three-dimensional (3D) adaptive grids. As promising results were obtained for turbulence-resolving techniques such as Direct Numerical Simulations and Large-eddy Simulation (LES), herein we explore whether similar advancements can be made with more practically oriented techniques for the numerical modelling of the atmosphere. As such, the present model uses Reynolds-averaged Navier-Stokes (RANS) techniques to parameterize the vertical mixing processes due to turbulence [START_REF] Reynolds | On the dynamical theory of incompressible viscous fluids and the determination of the criterion[END_REF], as is typical in weather and climate models.

1

The discussion of limited grid resolution is present in many studies of SCMs and GCMs. A prominent example is the nocturnal cumulus-cloud case (Wyant et al., 2007): whereas a high resolution mesh is required for capturing the processes at the cloud interface, a coarser resolution may be used for the time when the sun has risen and the cloud has been dissolved. More generally speaking, virtually all of the atmospheric dynamics that require a relatively high-resolution grid for their representation in numerical models are localized in both space and in time. The issue is made more difficult to tackle by the fact that their spatio-temporal localization is typically not known a priori (e.g. the height and strength of a future inversion layer). Therefore, the pre-tuned and static-type grids that most operational GCMs use (virtually all) are not flexible enough to capture all dynamical regimes accurately or efficiently. This also puts a large strain on the used closures for the sub-grid scale processes. In order to mitigate this challenge, GCMs that employ a so-called adaptive grid have been explored in the literature.

Here the grid resolution adaptively varies in both space and time, focussing the computational resources to where and when they are most necessary. Most notably, the innovative works of [START_REF] Jablonowski | Adaptive grids in weather and climate modeling[END_REF], [START_REF] Jablonowski | Block-structured adaptive meshes and reduced grids for atmospheric general circulation models[END_REF] and [START_REF] St-Cyr | A comparison of two shallow-water models with nonconforming adaptive grids[END_REF] report on the usage of both Cartesian and non-conforming three-dimensional adaptive grids and clearly demonstrate the potential of grid adaptivity for GCMs. Inspired by their works, we follow a 1D SCM approach and aim to add to their findings, using different grid-adaptative formulations and solver strategies. Since SCMs do not resolve large-scale atmospheric circulations, the analysis herein focusses on the representation of the Atmospheric Boundary Layer (ABL).

Over the years, the computational resources that are available to run computer models have increased considerably [START_REF] Schaller | Moore's law: past, present and future[END_REF]. This has facilitated GCMs to increase their models' spatial resolution, enabling to resolve the most demanding processes with increased grid resolutions. However, it is important to realize that the (spatial and temporal) fraction of the domain that benefits most from an increasing maximum resolution necessarily decreases as separation of the modelled spatial scales increases [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF]. This is because the physical processes that warrant a higher-resolution mesh are virtually never space filling. E.g. the formation phase of tropical cyclones is localized in both space and time and is characterized by internal dynamics that evolve during the formation process. By definition, with an increasing scale separation, only an adaptive-grid approach is able to reflect the effective (or so-called fractal) dimension of the physical system in the scaling of the computational costs [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF][START_REF] Van Hooft | Towards Adaptive Grids for Atmospheric Boundary Layer Simulations[END_REF]. This is an aspect where the present adaptive-grid approach differs from for example, a dynamic-grid approach [START_REF] Dunbar | A one-dimensional finite-element boundary-layer model with a vertical adaptive grid[END_REF], that employs a fixed number of grid cells that needs to be predefined by the user. This work employs a similar method for grid adaptation as presented in the work of [START_REF] Van Hooft | Towards Adaptive Grids for Atmospheric Boundary Layer Simulations[END_REF] on 3D-turbulence-resolving simulations of the ABL. As such, this work is also based on the adaptive-grid toolbox and built-in solvers provided by the 'Basilisk' code (http://basilisk.fr).

We test our model with the well established cases defined for the first two GABLS intercomparison projects for SCMs. As part of the Global Energy and Water cycle EXchanges (GEWEX) modelling and prediction panel, the GEWEX ABL Study (GABLS) was initiated in 2001 to improve understanding of the atmospheric boundary layer processes and their representation in models. Based on observations during field campaigns, a variety of model cases has been designed and studied using both LES and SCMs with a large set of models using traditional static-grid structures. An overview of the results and their interpretation for the first three intercomparison cases are presented in the work of [START_REF] Holtslag | Stable atmospheric boundary layers and diurnal cycles: challenges for weather and climate models[END_REF]. Here we will test the present adaptive-grid SCM based on the first two intercomparison cases, referred to as GABLS1 and GABLS2. These cases were designed to study the model representation of the stable boundary layer and the diurnal cycle, respectively. Their scenarios prescribe idealized atmospheric conditions and lack the complete set of physical processes and interactions encountered in reality. At this stage within our research, the authors consider this aspect to be an advantage, as the present SCM model does not have a complete set of parameterizations for all processes that are typically found in the operational models (see e.g. [START_REF] Slingo | The development and verification of a cloud prediction scheme for the ECMWF model[END_REF][START_REF] Grell | Fully coupled "online" chemistry within the WRF model[END_REF]). This paper is organized as follows, the present SCM is discussed in more detail in Sect. 2. Based on the results from a simplified flow problem, Sect. 3 starts with an analysis of the used numerical methods and the grid adaptation strategy. Model results for ABL-focussed cases that are based on the first two GABLS intercomparison scenarios are also presented in Sect. 3.

Finally, a discussion and conclusions are presented in Sect. 4.

Model Overview

As we focus on the merits of grid adaptivity in this study on SCMs and not on the state-of-the-art closures for the vertical transport phenomena, we have opted to employ simple and well-known descriptions for the turbulent transport processes. More specifically, the present model uses a stability-dependent, first-order, local, K-diffusivity closure as presented in the work of [START_REF] Louis | A short history of PBL parameterization at ECMWF[END_REF] and [START_REF] Holtslag | Local versus nonlocal boundary-layer diffusion in a global climate model[END_REF]. For the surface-flux parameterizations we again follow the formulations in the work of [START_REF] Holtslag | Local versus nonlocal boundary-layer diffusion in a global climate model[END_REF]. However, to improve the representation of mixing under stable conditions, an alteration is made to the formulation of the so-called stability-correction function under stably-stratified conditions. Based on the work of [START_REF] England | Stability functions based upon shear functions[END_REF], we use a so-called short-tail mixing function. The used closures for the turbulent transport are summarized next. The upward surface fluxes (F ) of the horizontal velocity components (u, v), the potential temperature (θ) and specific humidity q are evaluated as:

F u = -C M U 1 u 1 , (1a)
F v = -C M U 1 v 1 , (1b)
F θ = -C H U 1 (θ 1 -θ 0 ) , (1c)
F q = -C H U 1 (q 1 -q 0 ) , ( 1d
)
Where U is the wind-speed magnitude and indices 0 and 1 refer the to values at the surface and the first model level, respectively.

The surface transport coefficients are,

C M = C N f s,M (Ri b ), (2a)
C H = C N f s,H (Ri b ), (2b)
with Ri b the surface bulk Richardson number, that is defined as,

Ri b = g θ v,ref z 1 (θ v,1 -θ v,0 ) U 2 1 , (3)
where g is the acceleration due to gravity, θ v is the virtual potential temperature and θ v,ref is a reference temperature whose value is taken as a scenario-specific constant. Equation 3 assumes that θ v is related to the buoyancy (b) [START_REF] Boussinesq | Théorie de l'écoulement tourbillonnant et tumultueux des liquides dans les lits rectilignes à grande section[END_REF] via

g and θ v,ref according to b = g/θ v,ref (θ v -θ v,ref ) .
The virtual potential temperature is related to the potential temperature (θ) and specific humidity (q) according to,

θ v = θ 1 -1 - R v R d q , (4)
with R v /R d = 1.61 the ratio of the gas constants for water vapour and dry air [START_REF] Emauel | Atmospheric convection[END_REF][START_REF] Heus | Formulation of and numerical studies with the Dutch Atmospheric Large-Eddy Simulation (DALES)[END_REF]. The so-called neutral exchange coefficient (C N ) is calculated using,

C N = k 2 ln ((z 1 + z 0,M )/z 0,M ) 2 , (5)
with k = 0.4 the Von Karman constant, z 1 the height of lowest model level and z 0,M is the roughness length for momentum.

For the cases studied in this work, the roughness length for heat is assumed to be identical to z 0,M . The stability correction functions for the surface transport of momentum and heat (f s,M , f s,H ) are,

f s,M (Ri b ) =            0, Ri b ≥ 0.2, 1 -Ri b 0.2 2 , 0 ≤ Ri b < 0.2, 1 - 10Ri b 1+75C N √ ((z1+z 0,M )/z 0,M ) Ri b , Ri b < 0, (6a) f s,H (Ri b ) =      f s,M (Ri b ), Ri b ≥ 0, 1 - 15Ri b 1+75C N √ ((z1+z 0,M )/z 0,M ) Ri b , Ri b < 0, (6b)
which conclude the description of the surface fluxes. The vertical flux (w a ) of a dummy variable a due to turbulence within the boundary layer is based on a local diffusion scheme and is expressed as,

w a = -K ∂a ∂z , ( 7
)
where K is the so-called eddy diffusivity,

K = l 2 Sf (Ri). (8)
l represents an effective mixing length,

l = min (kz, l bl ) , (9)
with l bl is the Blackadar length scale, we use, l bl = 70m [START_REF] Holtslag | Local versus nonlocal boundary-layer diffusion in a global climate model[END_REF]. S is the local wind-shear magnitude,

S = ∂u ∂z 2 + ∂v ∂z 2 (10)
and f (Ri) is the stability correction function for the vertical flux,

f (Ri) =            0, Ri ≥ 0.2, 1 -Ri 0.2 2 , 0 ≤ Ri < 0.2, √ 1 -18Ri, Ri < 0, (11)
i.e. based on the gradient Richardson number,

Ri = g θ v,ref ∂θ v /∂z S 2 . ( 12
)
The authors of this work realize that there have been considerable advancements on the representation of mixing under unstable conditions in the past decades, e.g non-local mixing [START_REF] Holtslag | Local versus nonlocal boundary-layer diffusion in a global climate model[END_REF] and turbulent-kinetic-energy-based closures (see e.g., [START_REF] Mellor | Development of a turbulence closure model for geophysical fluid problems[END_REF][START_REF] Lenderink | An updated length-scale formulation for turbulent mixing in clear and cloudy boundary layers[END_REF]. Therefore, we would like to note that such schemes are compatible with the adaptive-grid approach and they could be readily employed to improve the physical descriptions in the present model. From an implementations' perspective, those schemes would not require any grid-adaptation specific considerations when using the Basilisk code.

For time integration; we recognize a reaction-diffusion-type equation describing the evolution of the horizontal wind components and scalar fields such as the virtual potential temperature and specific humidity (q). For a variable field s(z, t), we write,

∂s ∂t = ∂ ∂z (K ∂ ∂z s) + r. ( 13
)
Where r is a source term and K is the diffusion coefficient c.f. Eq. 8. Using a mixed implicit-explicit first-order-accurate time discretization for the diffusive term and an explicit time integration for the source term (r), with time step ∆t separating the solution s n and s n+1 , this can be written,

s n+1 -s n ∆t = ∂ ∂z (K n ∂ ∂z s n+1 ) + r n . ( 14
)
Rearranging the terms we write,

∂ ∂z (K n ∂ ∂z s n+1 ) - s n+1 ∆t = - s n ∆t -r n , (15)
to obtain a Poisson-Helmholtz equation for s n+1 , using the eddy diffusivity calculated from the solution s n (K n ). Eq. 15 is solved using a multigrid strategy, employing a finite-volume-type second-order-accurate spatial discretization (Popinet, 2017a, b). The source term r in Eq. 13 is defined using different formulations for the various scalar fields in our model. For θ and q, the source term r concerns the tendency in the lowest grid level due to the surface fluxes (F , see Eqs. 1, r Fs ) and the effect of large scale synoptic divergence (r w ) according to the vertical velocity w (i.e. prescribed for the GABLS2 case). We write for a dummy variable s,

r w,s = -w ∂s ∂z (16)
For the horizontal velocity components (u, v) the corresponding source terms (i.e. r Fs and r w ) are also taken into account and supplemented with the additional source term r ∇hP,f , that concerns the horizontal pressure-gradient-forcing vector (i.e.

-1 ρ ∇ h P , for air with a density ρ) and the Coriolis-force term according to the local Coriolis parameter f . For the horizontal velocity vector u = {u, v, 0} we write,

r ∇hP,f = -∇ h P ρ -f k × u , (17)
where '×' represents the cross product operator and k = {0, 0, 1} the unit vector in the vertical direction. In this work we adopt the commonly used strategy to introduce a velocity vector that known as the geostrophic wind (U geo ), according to,

U geo = k ρf × ∇ h P. (18)
The most prominent feature of the SCM presented in this work is that it adaptively coarsens and refines the grid resolution based on the evolution of the solution itself. As mentioned in the introduction, the associated grid-adaptation algorithm is the same as described in [START_REF] Van Hooft | Towards Adaptive Grids for Atmospheric Boundary Layer Simulations[END_REF]. Here we only briefly discuss the general concept.

Apart from the imperfect representation of the physical aspects of a system in numerical models, additional errors naturally arise due to the spatial and temporal discretization. In general, a finer resolution corresponds to a more accurate solution and a simulation result is considered to be 'converged' when the numerically obtained solution and the statistics of interest do not crucially depend on the chosen resolution. The aim of the grid-adaptation algorithm is to dynamically coarsen and refine the mesh so that the errors due to the spatial discretization remain within limited bounds and to be uniformly distributed in both space and time. For our adaptive approach this requires, (1) an algorithm that evaluates a local estimate of the discretization error in the representation of selected solution fields (χ a for a field 'a') and ( 2), a corresponding error threshold (ζ a ) that determines if a grid cell's resolution is either too coarse (i.e. χ a > ζ a ), too fine (i.e. χ a < 2ζ a /3), or just fine. Grid adaptation can then be carried out accordingly and the solution values of new grid cells can be found using interpolation techniques. A cell is refined when the estimated error for at least one selected solution field exceeds it's refinement criterion and a cell is coarsened when it is considered to be 'too fine' for all selected solution fields. The 'error estimator' (χ) is based on a so-called multi-resolution analysis that is formally linked to wavelet thresholding. The algorithm aims to estimate the magnitude of higher-order contributions in the spatial variability of the solution that are not captured by the solver's numerical schemes.

Consistent with the second-order spatial accuracy of the solver's numerical schemes [START_REF] Popinet | Multigrid Poisson-Helmholtz solvers[END_REF], the algorithm employs a second-order accurate wavelet-based error estimate. In practice, grid refinement will typically occur at the locations where the solution is highly 'curved', whereas those areas where the solution fields vary more 'linearly' in space are prone to coarsening.

The error threshold, or so-called refinement criterion ζ, is defined by the user. Noting that similar to the pre-tuning of the fixed-in-time grids as is common in most SCMs, the balance between accuracy and the required computational effort remains at the discretion of the model's user.

For the cases in this work that focus on the ABL (i.e. in Sect 3.2 and 3.3), the dynamics are governed by the wind (U = (u, v))

and the virtual potential temperature (θ v ), hence we base the refinement and coarsening of the grid on a second-order-accurate estimated error associated with the representation of these discretized fields. Based on trial and error, we set the corresponding refinement thresholds,

ζ u,v = 0.25 m/s, (19)
ζ θv = 0.5 K, (20)
for both of the horizontal wind components and virtual potential temperature, respectively. These values are the result of a choice by the authors that aims to strike an arbitrary balance between the accuracy of the solution and the computational effort required to run the model. Note that a similar (arbitrary) balance needs also to be found when static grids are employed. For a simple flow set-up, Sect. 3.1 presents an example convergence study to show the effects of using different refinement criteria on the accuracy of the obtained solutions.

Grid adaptation is carried out each time step. The tree-based anisotropic-grid structure in Basilisk facilitates a convenient basis for the multi-resolution analysis and the subsequent refinement and coarsening of cells at integer levels of refinement.

This entails that the spatial resolution can vary by factors of two [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF]. For the adaptive-grid runs presented in this paper, the time spent in the actual grid assessment and adaptation routines is less than than 5% of the total wall-clock time (see table 1).

Apart from the Ekman-spiral case in Sect. 3.1, the physical time step in the ABL-focussed cases is adaptively varied between 2 sec. and 15 sec. based on the convergence properties of the aforementioned iterative solver. Noting that these values are rather small compared to existing GCMs that often employ higher-order-accurate time-integration schemes. Additionally, the correlation of spatial and temporal scales warrants a smaller time step, since the present model employs a higher maximum vertical resolution compared to that of an operational GCM. The solver's second-order spatial accuracy is validated and the performance is accessed for a simple flow set-up in Sect. 3.1. For the exact details of the model set-ups for the cases presented in this paper, the reader is referred to the case-definition files (in legible formatting). Links are provided to their online locations in table 1.

The Laminar Ekman spiral and grid adaptation

Before we focus our attention on cases that concern the ABL, this section discusses the philosophy of the used grid adaptation strategy based on the analysis of a one-dimensional (1D) laminar Ekman-flow set-up. This simple and clean set-up enables to quantify numerical errors explicitly and test the solver's numerical schemes. The aim of this section is to show that the gridadaptation strategy and the accompanying refinement criteria provide a consistent and powerful framework for adaptive meshelement-size selection. Results are presented for both an equidistant-grid and the adaptive-grid approach. The case describes a neutrally-stratified fluid with a constant diffusivity for momentum (K) given by the kinematic viscosity ν and density ρ in a rotating frame of reference with respect to the Coriolis parameter f . A flow is forced by a horizontal pressure gradient (-∇ h P ) according to Eq. 18 using U geo = {U geo , 0}, over a no-slip bottom boundary (located at z bottom = 0). Assuming that the velocity components converge towards the geostrophic wind vector for z → ∞ and vanish at the bottom boundary, there exists an analytical, 1D, steady solution for the horizontal wind component profiles (u E (z), v E (z)), that is known as the Ekman spiral;

u E = U geo 1 -e -γz cos(γz) , (21)
v E = U geo e -γz sin(γz), (22)
with γ the so-called inverse Ekman depth, γ = f / (2ν). We choose numerical values for U geo , γ and f of unity in our set-up and present the results in a dimensionless framework. The solution is initialized according to the exact solution and we set boundary conditions based on Eqs. 21 and 22. Equation 13 is solved numerically for both u and v components, on a domain with height z top = 100γ -1 . The simulation is run until t end = 10f -1 , using a fixed time step ∆t = 0.01f -1 . The time step is chosen sufficiently small such that the numerical errors are dominated by the spatial discretization rather than by the timeintegration scheme. During the simulation run, discretization errors alter the numerical solution from it's exact, and analytically steady, initialization. For all runs, the diagnosed statistics regarding the numerical solutions that are presented in this section have become steady at t = t end .

The spatial-convergence properties for the equidistant-grid solver are studied by iteratively decreasing the used (equidistant) mesh-element sizes (∆) by factors of two and we monitor the increasing fidelity of the solution at t = t end between the runs.

Therefore, based on the analytical solution, a local error ( u,v ) of the numerically obtained solution (u n , v n ) within each grid cell is diagnosed and is defined here as:

a = a n -a E , ( 23
)
where a is a dummy variable for u and v, a E is the grid-cell-averaged value of the analytical solution (a E ) and a n the value of the numerical solution within the cell. Noting that a n also represents the grid-cell-averaged value in our finite-volume approach. Figure 1a shows the results for all runs and compares the used grid resolution (∆) with the error u,v . It appears that the observed range of -values is large and typically spans 10 orders of magnitude, with a lower bound defined by the 'machine precision' (i.e. ≈ 10 -15 for double-precision floating-point numbers). This wide range can be explained by the fact that the Ekman spiral is characterized by exponentially decreasing variation with height (see Eqs. 21,22) and hence the equidistant grid may be considered overly refined at large z. This illustrates that, for a given solver formulation, the error in the solution is not directly dictated by the mesh-element size, but also depends on the local shape of the numerical solution itself. This poses a challenge for the pre-tuning of meshes applied to GCMs, where a balance need to be found between accuracy and computational speed performance. The solution of a future model run is not known beforehand and hence the tuning of the grid typically relies heavily on experience, empiricism and a-priori knowledge. This motivates to apply the method of error estimation in the representation of a discretized solution field as described in [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF] and [START_REF] Van Hooft | Towards Adaptive Grids for Atmospheric Boundary Layer Simulations[END_REF].

For both velocity components, this estimated error (χ u,v ) is evaluated at the end of each simulation run for each grid cell and is plotted against the corresponding actual error ( u,v ) in Fig. 1b. It seems that for this virtually steady case, there is a clear correlation between the diagnosed (instantaneous) χ-values and the -values that have accumulated over the simulation run time. Even though the correlation is not perfect, it provides a convenient and consistent framework for a grid adaptation algorithm. As such, a second convergence test for this case is performed using a variable-resolution grid within the domain. The mesh is based on the aforementioned adaptive-grid approach. For these runs, we iteratively decrease the so-called refinement criterion (ζ u,v ) by factors of two between the runs and monitor the increasing fidelity of the numerically obtained solution for all runs. The refinement criterion presets a threshold value (ζ) for the estimated error χ that defines when a cell should be refined

(χ > ζ)
or alternatively, when it may be coarsened (χ < 2ζ/3). Figure 2a presents the results and compares the used local grid resolution against u,v for the various (colour-coded) runs. It appears that for all separate runs, the algorithm employed a variable resolution mesh and that this has resulted in a smaller range of the local error in the solution ( ), as compared to the equidistant-grid cases. The local error in the solution is also compared against the wavelet-based estimated error in the representation of the solution fields in Fig. 2b. Compared to the results from the equidistant-grid approach as presented in Fig. 1b), the spread of the χ and values is relatively small for the separate runs when the adaptive-grid approach is used. The most prominent reason for the finite spread is that the error ( ) was diagnosed after 1000 time steps. This facilitated errors in the solution that arise in the solution at a specific location (with a large χ-value) to 'diffuse' over time towards regions where the solution remains to be characterized by a small χ-value (not shown). Also, since u and v are coupled (due to the background rotation), local errors that arise in the solution for u 'pollute' the v-component solution, and vice versa. Furthermore, a spread is expected because the tree-grid structure only allows the resolution to vary by factors of two [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF].

Finally, the global convergence characteristics and the speed performance of the two approaches are studied. The global error (η) in the numerically obtained solution is evaluated as,

η = ztop z bottom ( u + v ) dz, (24)
In order to facilitate a comparison between the methods, we diagnose the number of used grid cells (N ) for the adaptive-grid run. Figure 3a shows that for both approaches the error scales inversely proportional to the used number of grid cells to the second power (i.e. second-order spatial accuracy in 1D). The adaptive grid results are more accurate than the results from the fixed-grid approach when employing the same number of grid cells. Figure 3b shows that for both approaches the required effort (i.e. measured here in wall-clock time) scales linearly with the number of grid cells, except for the runs that require less than one-tenth of a second to perform. The plots reveals that per grid cell there is computational overhead for the adaptive-grid approach. These results show that the used numerical solver is well behaved.

The following sections are devoted to testing the adaptive-grid approach in a more applied SCM scenario, where the turbulent transport closures are applied (see Sect. 2) and the set-up is unsteady. Here, the quality of the adaptive-grid solution has to be assessed by comparing against reference results from other SCMs, large-eddy simulations and the present model running in equidistant-grid mode.

GABLS1

The first GABLS intercomparison case focusses on the representation of a stable boundary layer. It's scenario was inspired by the LES study of an ABL over the Arctic sea by [START_REF] Kosović | A large eddy simulation study of a quasi-steady, stably stratified atmospheric boundary layer[END_REF]. The results from the participating SCMs are summarized and discussed in [START_REF] Cuxart | Singlecolumn model intercomparison for a stably stratified atmospheric boundary layer[END_REF], for the LES intercomparison study, the reader is referred to the work of [START_REF] Beare | An intercomparison of large-eddy simulations of the stable boundary layer[END_REF]. The case prescribes the initial profiles for wind and temperature, a constant forcing for momentum corresponding to a geostrophic wind vector, U geo = {8, 0}m/s and Coriolis parameter f =1.39×10 -4 s -1 . Furthermore, a fixed surface-cooling rate of 0.25 K/hour is applied and θ v,ref = 263.5K. The model is run with a maximum resolution of 6.25 meter and a domain height of 400 meters. The maximum resolution corresponds to 6 levels of tree-grid refinement, where each possible coarser level corresponds to a factor of two increase in grid size.

Due to the idealizations in the case set-up with respect to the reality of the field observations, the model results were not compared against the experimental data [START_REF] Cuxart | Singlecolumn model intercomparison for a stably stratified atmospheric boundary layer[END_REF]. However, for the SCMs, a reference was found in the highfidelity LES results that tended to agree well between the various models. The LES results therefore serve as a benchmark for the results obtained with the present model. This facilitates a straightforward testing of the formulations and implementations of the used physical closures, before we continue our analysis towards the full diurnal cycle. Inspired by the analysis of [START_REF] Cuxart | Singlecolumn model intercomparison for a stably stratified atmospheric boundary layer[END_REF] and their figure 3, we compare our SCM results with the 6.25 meter-resolution LES ensemble results. We focus on the profiles for the wind components and potential temperature averaged over the ninth hour of the simulation in Fig. 4.

We observe that the present SCM is in good agreement with the LES results and is able to capture the vertical structure of the ABL, including the low-level jet. The differences are only minor compared to the variations in the results presented in the aforementioned GABLS1 SCM reference paper.

Note that in general, results are of course sensitive to the closure chosen to parameterize the turbulent transport, in our case given by Eqs.6 and 11. In order to separate between the numerical effects of using grid adaptivity and the chosen physical closures, we define an additional reference case in which we run an equidistant-grid SCM. This model run employs a fixed 6.25 meter resolution (i.e with 64 cells), but otherwise identical closures and numerical formulations. I.e. we have switched-off the grid adaptivity and maintain the maximum resolution throughout the domain. We can observe that results between both SCMs are in good agreement but that minor deviations are present. These discrepancies are on the order of magnitude of the refinement criteria and can be reduced by choosing more stringent values, that would result in using more grid cells. The evolution of the adaptive-grid structure is shown in Fig. 5 a. We see that a relatively high resolution is employed near the surface, i.e. in the logarithmic layer. Remarkably, without any a priori knowledge, the grid is refined at a height of 150 m < z < 200 m as the so-called low-level jet develops, whereas the grid has remained coarse above the boundary layer where the grid resolution was reduced to be as coarse as 100 meters. From Fig. 5 b we learn that the number of grid cells varied between 11 and 24 over the course of the simulation run.

GABLS2

The second GABLS model intercomparison case was designed to study the model representation of the ABL over the course of two consecutive diurnal cycles. The case is set-up after the observations that were collected on the 23rd and 24th of October, 1999 during the CASES-99 field experiment in Leon, Kansas, USA [START_REF] Poulos | CASES-99: A comprehensive investigation of the stable nocturnal boundary layer[END_REF]. The case prescribes idealized forcings for two consecutive days that were characterized by a strong diurnal cycle pattern. During these days, the ABL was relatively dry, there were few clouds and θ v,ref = 283.15K. The details of the case are described in the work of [START_REF] Svensson | Evaluation of the diurnal cycle in the atmospheric boundary layer over land as represented by a variety of single-column models: the second GABLS experiment[END_REF] that was dedicated to the evaluation of the SCM results for the GABLS2 intercomparison. Compared to the original case prescriptions, we choose a slightly higher domain size of z top = 4096 meters (compared to 4000 m), so that a maximum resolution of 8 meters corresponds to 9 levels of refinement.

In this section we place our model output in the context of the results presented in the work of [START_REF] Svensson | Evaluation of the diurnal cycle in the atmospheric boundary layer over land as represented by a variety of single-column models: the second GABLS experiment[END_REF], that, apart from the SCM results, also includes the results from the LES by [START_REF] Kumar | Impact of surface flux formulations and geostrophic forcing on large-eddy simulations of diurnal atmospheric boundary layer flow[END_REF]. To obtain their data we have used the so-called 'data digitizer' of [START_REF] Rohatgi | WebPlotDigitizer[END_REF]. Inspired by the analysis of [START_REF] Svensson | Evaluation of the diurnal cycle in the atmospheric boundary layer over land as represented by a variety of single-column models: the second GABLS experiment[END_REF] and their figures 10 and 11, we compare our results for the wind-speed magnitude (U = u ) and virtual potential temperature profiles at 14:00 local time on the 23rd of October in Fig. 6 a and b, respectively. Here we see that the results obtained with the present SCM fall within the range of the results as were found with the selected models that participated in the original intercomparison. These models also employed a first-order-style turbulence closure and have a lowest model-level height of less than 5 meters. The present modelled virtual potential temperature (θ v ) shows a slight negative vertical gradient in the well mixed layer. This is a feature related to the usage of the local K-diffusion description for the turbulent transport (see Sect. 2 and the work of [START_REF] Holtslag | Local versus nonlocal boundary-layer diffusion in a global climate model[END_REF]). Figure 6 c presents a timeseries of the 10-meter wind speed (U 10m ) during the 23-rd of October. Again the present model results compare well with the others. Next, in order to validate the grid-adaptivity independently from the used closures, we present the hourly evolution of the wind speed on the 24-th of October against the results obtained with adaptivity switched off, using 512 equally-spaced grid points in Fig. 7. A nearly identical evolution of the wind speed profiles is observed and even the small-scale features in the inversion layer (i.e. z ≈ 800 m) are present in the adaptive-grid-model calculations.

The corresponding evolution of the adaptive-grid structure is presented in Fig. 8, where the colours in the resolution plot appear to sketch a 'Stullian' image, showing a prototypical diurnal evolution of the ABL (see figure 1.7 in the book of [START_REF] Stull | An introduction to boundary layer meteorology[END_REF]. Apparently, the grid-adaptation algorithm has identified (!) the 'surface layer' within the convective boundary layer, the stable boundary layer, the entrainment zone and the inversion layer as the dynamic regions that require a high-resolution mesh.

Conversely, the well-mixed layer within the CBL, the residual layer and the free-troposphere are evaluated on a coarser mesh.

The total number of grid cells varied between 24 and 44.

experiments were set up and performed by AvH. The writing was led by AvH and organized in an iterative procedure with BvdW.

16 The wall-clock times are evaluated using a single core (processor model: Intel i7-6700 HQ). for 20 runs using the adaptive-grid approach with different refinement criteria (see colour bar). The left-hand side plot (a) Shows that the diagnosed errors for each run plotted against the used mesh-element size (∆). The inset (using the same axis scales) shows the results for a single run. The right-hand side plot (b) shows the correlation between the wavelet-based estimated error (χ) and the corresponding diagnosed error in the numerically obtained solution ( )). The inset (using the same axis scales) shows the results for a single run, and reveals a relatively small spread in both and χ values compared to the equidistant-grid results presented in figure 1b. 2006) (grey-shaded areas) and the present SCM, employing an equidistant and static grid with a 6.25 meter resolution (dashed lines). For z > 250m, the profiles have remained as they were initialized.

Figure 1 .

 1 Figure 1. The locally evaluated error in the numerical solutions for u and v at t = t end , based on the analytical solution ( u,v , see Eq. 23) for 10 runs using different equidistant mesh-element sizes. The left-hand side plot (a) Shows that the diagnosed errors for each run plotted against the used mesh-element size (∆) times the inverse Ekman depth (γ, see text). The right-hand side plot (b) shows, with the same colour coding as in the left-hand side plot (a), the correlation between the wavelet-based estimated error (χ) and the corresponding diagnosed error in the numerically obtained solution ( ). The inset (using the same axis scales) shows the results for a single run, and reveal a spread of several orders of magnitude in both and χ values.

Figure 2 .

 2 Figure 2. The locally evaluated error in the numerical solutions for u and v at t = t end , based on the analytical solution ( u,v , see Eq. 23)

Figure 3 .

 3 Figure 3. The scaling characteristics for the laminar-Ekman-spiral case. (a) Presents the error convergence for the equidistant-grid and adaptive-grid approach. The errors (η) follow the slope of the blue dashed line that indicates the second-order accuracy of the methods. The wall-clock time for the different runs is presented in (b), showing that for both of the aforementioned approaches, the required effort scales linearly with the number of grid cells.

Figure 4 .

 4 Figure 4. Time averaged profiles over the ninth hour of the run according to the GABLS1 intercomparison scenario. For (a) the horizontal wind components and (b) the potential temperature. Results are obtained with the present adaptive-grid SCM (coloured lines), the LES models ensemble (i.e mean ± σ) from Beare et al. (2006) (grey-shaded areas) and the present SCM, employing an equidistant and static grid

Figure 5 .

 5 Figure 5. Evolution of (a) the vertical spatial-resolution distribution and (b) the total number of grid cells, for the GABLS1 intercomparison case.

Table 1 .

 1 The exact formulation of the methods are described at the online locations of the definition files for the different cases presented in

	this manuscript.
	Section Case	Grid	URL: http://www.basilisk.fr...	Number of time steps Wall-clock time
	3.1	Ekman spiral Adaptive	/sandbox/Antoonvh/ekman.c	1000 (×20 runs)	≈ 19 sec.
	"	"	Fixed & Equidistant /sandbox/Antoonvh/ekmanfg.c	1000 (×10 runs)	≈ 18 sec.
	3.2	GABLS1	Adaptive	/sandbox/Antoonvh/GABLS1.c	16204	≈ 1.4 sec.
	"	"	Fixed & Equidistant /sandbox/Antoonvh/GABLS1fg.c 16324	≈ 0.9 sec.
	3.3	GABLS2	Adaptive	/sandbox/Antoonvh/GABLS2.c	24262	≈ 9 sec.
	"	"	Fixed & Equidistant /sandbox/Antoonvh/GABLS2fg.c 33993	≈ 22 sec.

the GNU/Linux project and/or the GNU compiler collection (see www.gnu.org).

running the scripts. Finally, a snapshot of the used code, as it was used in this the work, is made available via ZENODO, with doi link: https://doi.org/10.5281/zenodo.1203631.

Acknowledgements. The authors gratefully acknowledge the funding by the ERC Consolidator grant (648666). The LES ensemble results used for the GABLS1 intercomparison are kindly made available by Bob Beare; online via: http://gabls.metoffice.com/ and we thank three anonymous reviewers for their comments on the manuscript. Furthermore, we acknowledge that we are indebted to those who contributed to

In this work we have presented a one-dimensional (1D) single-column model (SCM) that employs a mesh whose resolution is varied adaptively based on the evolution of the numerically obtained solution. This is an attractive feature because it is a prerequisite to enable the computational effort required for the evaluation of numerical solution to scale with the complexity of the studied physical system. The adaptation algorithm based on limiting discretization errors appears to function very well for a wide variety of geophysical applications: e.g. 3D atmospheric turbulence-resolving models [START_REF] Van Hooft | Towards Adaptive Grids for Atmospheric Boundary Layer Simulations[END_REF], tsunami and ocean-wave modelling [START_REF] Popinet | Quadtree-adaptive tsunami modelling[END_REF][START_REF] Beetham | Wave transformation and shoreline water level on Funafuti Atoll, Tuvalu[END_REF][START_REF] Marivela-Colmenarejo | Numerical Perspective on Tsunami Hazards and Their Mitigation by Coastal Vegetation[END_REF], hydrology [START_REF] Kirstetter | Modeling rain-driven overland flow: Empirical versus analytical friction terms in the shallow water approximation[END_REF], two-phase micro physics [START_REF] Howland | It's Harder to Splash on Soft Solids[END_REF], flow of granular media [START_REF] Zhou | Experiments on, and discrete and continuum simulations of, the discharge of granular media from silos with a lateral orifice[END_REF] and shock-wave formation [START_REF] Eggers | Spatial structure of shock formation[END_REF]. For these studies on highly dynamical systems, the adaptive-grid approach is chosen because it offers a more computationally efficient approach as compared to the usage of static grids.

The present work does not include an in-depth assessment and discussion on the performance of the presented methods in relation to the computational speed. Even though this is an important motivation for the application of the adaptive-grid strategy to GCMs, the authors argue that a SCM is not suitable for speed-performance testing: the speed of single-column calculations is virtually never a critical issue. Only in 3D mode, when SCMs are 'stitched together' to enable the resolving of global circulations, the model's computational efficiency becomes an issue. Furthermore, the performance of a SCM that employs a few tens of cells is not a good indicator for the performance of a GCM that can employ billions of grid cells. For the latter, parallel computation overhead and the so-called memory bottle neck are important aspects. In contrast, for the SCM case, the complete instruction set and solution data can typically be loaded onto the cache memory of a single CPU's core. Nevertheless, for the readers' reference, the required run times for the different SCM set-ups presented herein are listed in table 1, and figure 3b also presents quantitative results on this topic and shows that the adaptive-grid solver is well behaved.

Following the turbulence resolving study of [START_REF] Van Hooft | Towards Adaptive Grids for Atmospheric Boundary Layer Simulations[END_REF], the results presented herein are a proof-of-concept for future 3D modelling, using RANS techniques. The authors of this work realize that the present SCM is a far cry from a complete global model and that more research and development is required before the method can be compared on a global-circulation scale. As shown by e.g. [START_REF] Jablonowski | Adaptive grids in weather and climate modeling[END_REF], a 3D adaptive grid also allows a variable grid resolution in the horizontal directions. This further enables the computational resources to focus on the most challenging atmospheric processes where there is a temporal and spatial variation in the horizontal-resolution requirement of the grid. Examples include the contrasting dynamics between relatively calm centres of high-pressure circulations and those characterizing stormy low-pressure cells. Also, in the case of a sea breeze event [START_REF] Arrillaga | A characterisation of sea-breeze events in the eastern Cantabrian coast (Spain) from observational data and WRF simulations[END_REF], it would be beneficial to temporarily increase the horizontal resolution near the land-sea interface. As such, we encourage the usage of this technique for those meteorologically challenging scenarios.

Code and data availability. Basilisk is a freely available (GPLv3), multi-purpose tool to solve partial differential equations and has it's own website: http://www.basilisk.fr. The code contains solvers for Saint-Venant problems, the Navier-Stokes equations, electro hydrodynamics and more, see http://basilisk.fr/src/README. A selection of examples can be viewed online: http://www.basilisk.fr/src/examples. The website also provides general information including; installation instructions and a tutorial. Furthermore, for the work herein, interested Author contributions. All Authors contributed to the content of this manuscript and it should be viewed as a fruit of the many discussions we have had. Furthermore, SP wrote the Basilisk code, he also designed and implemented the grid adaptation algorithm. The numerical Competing interests. The authors declare that there are no competing interests. full diurnal cycles, corresponding to the 23rd and 24th day of October, 1999 (ranging from the labels 1:00:00 to 3:00:00 on the x-axis).


H Ferrage

D Taubald

B Schild

Velde

A Bauer

B Lanson

K Emmerich

E Ferrage

H Taubald

D Schild

B Velde

THE FATE OF SMECTITE IN KOH SOLUTIONS

whether they are published or not. The documents may come

Introduction

In 1997 at the end of a workshop dedicated to Victor Drits by the French Clay Group Bernard Kübler proposed that it would be a great help if the speakers would give their definitions of vermiculite, high-charge-smectite and illite at the beginning of their talks. As a matter of fact a look through literature reveals that the terms used for expandable layers in the smectite to illite reaction are just operationally defined and as such they do not necessarily coincide with definitions from the Nomenclature committee of the AIPEA. The aim of our study was to investigate the fate of smectite during a possible smectite to illite transformation. The question is then, is smectite just a passive partner in the reaction giving the ingredients for the new phase illite or is illite perhaps smectite with a higher charge and smectite collapse occurs only when potassium is supplied to the high-charge interlayer sites? Chemical compositions of diagenetically altered shales, sandstones or bentonites studied by [START_REF] Srodon | Chemistry of illite-smectite inferred from TEM measurements of fundamental particles[END_REF][START_REF] Srodon | Chemistry of Illite Smectite and End-Member Illite[END_REF][START_REF] Nadeau | Composition of Some Smectites and Diagenetic Illitic Clays and Implications for Their Origin[END_REF][START_REF] Awwiller | Illite/smectite formation and potassium mass transfer during burial diagenesis of mudrocks: a study from the Texas Gulf Coast Paleocene-Eocene[END_REF][START_REF] Cuadros | Characterization of mixed-layer illite-smectite from bentonites using microscopic, chemical, and X-ray methods: Constraints on the smectiteto-illite transformation mechanism[END_REF] indicate the formation of high-charge expandable layers together with the illitization process. The experimental illitization of smectite indicates the proportion of highcharge layers increases followed by the potassium fixation and the subsequent illitization [START_REF] Whitney | Experimental investigation of the smectite to illite reaction: Dual reaction mechanisms and oxygen-isotope systematics[END_REF]. These high-charge but exchangeable layers are comparable to the vermiculitic component described by [START_REF] Shutov | On the mechanism of a postsedimentary transformation of montmorillonite into hydromica[END_REF]Drits et al., 1997a. As a result, a randomly interstratified illite/smectite (R0 I/S) mixed-layer structure (MLS) has three components: one fully expandable (smectite), one partially expandable (vermiculite) and one non expandable (illite). This means that the smectite-to-illite transformation must involve the vermiculitization of the smectite layers through the formation of high-charge, tetrahedrallysubstituted layers from the original clay [START_REF] Drits | Mixed-layer minerals: Diffraction methods and structural features[END_REF]. Meunier et al., 2000 described this process as the vermiculitization of smectite interfaces. In their model the low-charge montmorillonite is transformed into a high-charge vermiculite which remains partly expandable (one water layer or one ethylene glycol layer at 80 % relative humidity) if the interlayer is not K + -or NH 4 + -saturated. A vermiculite interlayer leads to the formation of illite when K or NH 4 + is supplied. But how are smectite particles transformed during their transition to R0 I/S? Two processes may control the reaction by solid state transformation (Dunoyer de [START_REF] Dunoyer De Segonzac | The transformation of clay minerals during diagenesis and low-grade metamorphism: A review[END_REF][START_REF] Hower | Mechanism of burial metamorphism of argillaceous sediment: 1. Mineralogical and chemical evidence[END_REF] or dissolution/crystallisation (Baronnet, 1997, Altaner and[START_REF] Altaner | Comparison of structural models of mixed-layer illite/smectite and reaction mechanisms of smectite illitization[END_REF]. Many studies have reported evidence supporting one or the other reaction mechanism.

In the literature two experimental pathways are commonly used to obtain illite from a smectite precursor. One is to use KCl and high temperatures (> 200°C) and pressures and the other is to use KOH at low temperatures (< 100°C). Hydrothermal synthesis experiments have formed R0 I-S from smectite at elevated temperatures (e.g. [START_REF] Robertson | Smectite to illite conversion rates: Effects of solution chemisrty[END_REF]Lahann, 1981, Whitney and[START_REF] Whitney | Experimental investigation of the smectite to illite reaction: Dual reaction mechanisms and oxygen-isotope systematics[END_REF]. Reactions in KOH solutions lead to the decrease of the coherent scattering domain size of the smectite followed by the illitization of the smectite and finally the formation of pure illite (e. g. [START_REF] Eberl | Potassium Fixation in Smectite by Wetting and Drying[END_REF][START_REF] Eberl | Synthesis of illite-smectite from smectite at Earth surface temperatures and high pH[END_REF][START_REF] Bauer | Smectite transformation in high molar KOH solutions[END_REF]. In the present study we focus on the KOH system. The aim of the present manuscript is to follow the fate of a fully expandable and low-charge smectite before it becomes illite in KOH solutions with different experimental techniques.

Materials and Methods

Starting Clay: The starting clay mineral for this experimental study was purified homoionic montmorillonite obtained from the SAz-1 source clay material (Clay Minerals Society, Purdue University, http://www.clays.org/sourceclays/SourceClays.html). In a first step the clay was conditioned to remove soluble salts, sparingly soluble minerals and trace cation impurities after the procedure proposed by [START_REF] Bauer | Smectite transformation in high molar KOH solutions[END_REF]. The purified smectite was converted into its sodium form [START_REF] Moore | X-Ray Diffraction and the Identification and Analysis of Clay Minerals[END_REF] and the < 1 µm fraction of smectite was separated by sedimentation techniques [START_REF] Day | Particle Fractionation and particle size analysis[END_REF]. XRD analysis for conditioned clay revealed that montmorillonite is present at 100 %. Finally the clay separate was converted into its potassium form.

Starting solution:

The solutions were prepared from analytical grade chemicals. A Milli-Q Reagent Water System from Millipore Corp. provided deionized water (DI) with a resistivity of > 18 MΩ cm -1 .

Experimental protocol: Smectite samples weighing 2 g were added to 160 ml of 1 M KOH aqueous solution. Reactions were conducted in HDPE Nalgene bottles kept in ovens at 80°C

with reaction times lasting up to 150 days. The bottles were shaken twice a week. At the end of an experimental run the bottles were removed, quenched in cold water (25°C) and opened.

The supernatant was filtered through a 0.1 µm Nucleopore membrane into cleaned polypropylene bottles and analysed. The solid material was then washed carefully five times with 150 ml of DI water. After washing the material, it was five times resuspended for 1 minute in a ultrasonic bath in 150 ml of 1 M SrCl 2 solution. Saturation with Sr 2+ was used to ensure the presence of two water layers in the expandable layers under XRD data collection conditions (see below). To promote complete cation exchange, the suspensions were stored after the ultrasonic treatment for 4 h at 50°C. After this procedure the sample was washed five times with 150 ml of DI water. After washing the samples one part of the solid was used to prepare the XRD slide while the other rest was converted into the sodium form and freeze dried for layer charge measurements, oxygen isotope and XPS measurements as well as for thermal analysis.

X-ray diffraction:

The evolution of clay mineral structure and the formation of secondary phases were monitored by XRD. Oriented slides were prepared by pipetting a slurry of Srsaturated samples on a glass slide and drying it at room temperature for a few hours to obtain an air-dried preparation. XRD scans were recorded on a Bruker D5000 diffractometer (Cu radiation at 40 kV and 40 mA) equipped with an Ansyco rh-plus 2250 humidity control device coupled to an Anton Paar TTK450 chamber. The measuring range was from 2° to 50°2θ with a step size of 0.04° and 6 second counting time. The divergence slit, the two Soller slits, the antiscatter slit and the resolution slit were 0.5°, 2.3°, 2.3°, 0.5° and 0.06°, respectively. All XRD analyses were conducted at a relative humidity of 40% after an homogenisation period of about 15 min before the measurement.

The 001 reflection profile of SAz-1 smectite was modelled using the algorithms developed initially by [START_REF] Sakharov | Mixed-layer kaolinite-montmorillonite: A comparison of observed and calculated diffraction patterns[END_REF][START_REF] Drits | X-Ray structure analysis of mixed-layer minerals[END_REF], Drits et al., 1997a, Sakharov et al., 1999a MgAlO 4 (99.985% MgAlO 4 powder, Alfa Aesar) samples were purified with ethanol and afterwards sputtered with Ar ions for 10 to 80 sec. Based on the measurements of the ultrapure solids, the PHI Multipak sensitivity factors were refined for the Al K a X-ray excitation.

The detection limit of the Si, Al and Mg analyses is 0.5 Mol%.

Thermal Analysis:

The measurements were performed on a STA 449 C Jupiter from NETZSCH-Gerätebau GmbH with a TG/DSC sample holder. The STA is connected with a quadruple mass spectrometer 403 C Aëolos from IPI / InProcess Instruments / NETZSCH-Gerätebau GmbH. All samples were allowed to equilibrate at a relative humidity of 53% in a dessicator above saturated Mg(NO 3 ) 2 solution for at least 48 h. The conditions for the STA measurement are reported in Table 1.

Infra-red spectroscopy:

The IR spectra were recorded with a Perkin-Elmer 1760 spectrometer over the range 4000-400 cm -1 (optical resolution, 2 cm -1 ). The pellets were prepared by diluting 1 mg of sample in 100 mg of KBr.

AFM measurements: Both a raw sample of SAz-1 and samples of SAz-1 recovered from experiments were observed by AFM. A dilute suspension of each SAz-1 Na-saturated sample was gently sonicated and centrifuged onto a crystal of highly ordered pyrolytic graphite (HOPG). The excess aqueous solution was removed and the sample dried in an oven at 40°C.

A Veeco Instruments DIMENSION 3100 AFM with Nanoscope IV controller operating in tapping mode was used for imaging the smectite particles under ambient conditions in air (after drying). High aspect ratio Si single-crystal AFM tips were used. The nominal oscillation resonance frequency of the cantilever was 283 kHz. The maximum scan rate was 5 Hz. The smectite particle morphology and dimensions were quantitatively evaluated using the Image WSXM 2.0 software (NANOTEC). Using AFM, we analyzed smectite particles laying flat on the HOPG substrate. Aggregates of randomly oriented smectites were disregarded in the measurements.

Scanning electron microscopy: Scanning electron microscopy (SEM) was used to observe material of which there was too little to analyse by XRD. The SEM photographs were obtained using a JEOL GSM-6100 instrument with a KEVEX energy dispersive spectrometer used for semi-quantitative analysis (peak heights ratios).

Oxygen Isotope Analyses of clays:

The oxygen isotope composition of the smectite samples was measured at the University of Tübingen, using a conventional extraction of oxygen and BrF 5 reagent, according to a method adapted from [START_REF] Clayton | The use of bromine pentafluoride in the extraction of oxygen from oxides and silicates for isotopic analysis[END_REF]. About 7 mg of sample were used and the extracted oxygen was converted to CO 2 . The sample is wrapped in aluminium foil and dried for 3 days in an exsiccator with P 2 O 5 to remove non-structural water (Savin & Epstein, 1970). Then the samples are immediately (within seconds) loaded into the Ni-reaction vessel in a dry N 2 stream and pumped for at least 2 hours at temperatures below 180˚C before BrF 5 was added. Reaction was performed at 550ºC for 16-18 h. Isotope measurements were made using a Finnigan MAT isotope ratio mass spectrometer (Finnigan MAT 252). Oxygen isotope compositions are given in the standard δ-notation, expressed relative to VSMOW in permil (‰). The precision of δ 18 O values was better than 0.2‰ when compared to the reference δ 18 O values for NBS-28 of 9.64‰.

The 1 M KOH fluid was analysed for total oxygen isotope composition (H 2 O + KOH) using a TC-EA preparation device connected online with a Finnigan MAT Delta + XL continuous flow mass-spectrometer. 0.5 µl liquid was filled in a silver container and combusted at 1450°

C. The analytical precision is better than ± 0.3‰.

Solution compositions: Concentrations of dissolved Si, Al, K, Fe and Mg were determined for all samples. The Si and Al were analysed using a Plasma 400 ICP-AES (Perkin Elmer).

All measurements were done in triplicate.

Results

XRD:

For Sr-saturated SAz-1 specimen the position of the 001 reflection is ~15.3 Å (Fig. 1).

With increasing reaction time, the profile of this reflection is slightly broadened and becomes asymmetrical on the high-angle side. Profile broadening may be related to a decrease of the coherent scattering domain size but the asymmetry on the high-angle side of the 001 reflection is more likely attributed to the presence of layers with lower layer thickness. The increased proportion of these layers is best revealed by XRD profile modelling. Structure models and structural parameters used to fit the 001 reflection are reported in Table 2. In agreement with structure models determined by [START_REF] Ferrage | Etude expérimentale de l'hydratation des smectites par simulation des raies 00L de diffraction des rayons X. Implications pour l'étude d'une perturbation thermique sur la minéralogie de l'argilite du site Meuse-Haute Marne[END_REF] on the <1 µm fraction of Sr-saturated SAz-1, two MLSs were used to reproduce the 001 reflection. The first one is mainly constituted by 2W layers with minor 1W layers (S1) whereas the second one is more heterogeneous and incorporates the three layer types (S2).

XRD measurements at 40 % RH of the Sr-saturated samples reveal some changes in the expandability of the smectite. With increasing alteration time the compositions of the two MLSs remain constant except for the last two samples in which an increase of 1W layers is observed in both contributions. A modification in expandability is observed, mainly induced by the increase of the S2:S1 ratio for reaction times up to 45 days (Table 2). As a result the number of 1W layers increases whereas that of 0W layers seems to be about constant with reaction time. The layer thickness of the different layer types were found constant whereas CSDS of the crystals decreases from 5.4 to 4.3 after 45 days and remains then constant till the end of the experiment (Table 2). Some SrCl 2 and SrCO 3 peaks were found as impurities resulting from the Sr 2+ -saturation procedure. In addition, low intensity sharp peaks, likely attributed to zeolite-type minerals (phillipsite), were observed for reacted samples.

Observation of the ( 060) peak region gives information on the a-b portion of the phyllosilicate structures. This region was monitored as a function of reaction time (Fig. 2). The appearance of a small peak at 1.54 Å beside the initial peak at 1.49 Å indicates the formation of trioctahedral coherently diffracting domains [START_REF] Brindley | Crystal Structures of Clay Minerals and their X-ray Identification[END_REF]. A similar observation was made by [START_REF] Beaufort | An experimental alteration of montmorillonite to a di plus trioctahedral smectite assemblage at 100 and 200 degrees C[END_REF] under neutral pH conditions at higher temperatures (100° & 200°C). More importantly, there is a clear shift in the position of the (060) peak of the major phase towards higher spacings. Is indicates a change in the internal structure of the aluminous smectite. Such a shift is indicative of an increase in layer charge, either for beidellite or montmorillonite substitutions in the lattice [START_REF] Velde | Clays Minerals: A physico-chemical explanation of their occurence[END_REF].

Solution chemistry:

The changes in solution composition with time were monitored during the reaction. Due to the high KOH content, no significant variation of pH has been observed.

After a first strong increase of the Al and Si concentration in solution the concentration increase flattens out to remain constant within the error of the measurement at ~ 950 mg/L for Si and ~ 250 mg/L for Al (Fig. 3). With increasing reaction time the Al/Si ratio in solution decreases from 0.39 to 0.25. The solution concentrations of Mg and Fe were about 0.06 mg/L and 0.3 mg/L, respectively, throughout the experiment.

The pH and aqueous species, as well as the saturation index of possible precipitates, were calculated at in-situ conditions by using the EQ3/6 software package [START_REF] Wolery | The Eq3/6 Software Package for Geochemical Modeling -Current Status[END_REF] and the SUPCRT92 database [START_REF] Johnson | Supcrt92 -a Software Package for Calculating the Standard Molal Thermodynamic Properties of Minerals, Gases, Aqueous Species, and Reactions from 1-Bar to 5000-Bar and 0-Degrees-C to 1000-Degrees[END_REF]. The in-situ pH was determined by calculating the charge balance on H + . We checked that the presence of small amounts of other cations does not change the pH calculations, given the high potassium concentrations.

Layer Charge:

We observed a continuous evolution of the layer charge with increasing reaction time. We measured a mean layer charge of 0.32 eq/(Si,Al) 4 O 10 for the initial SAz-1 which is in agreement with the value reported by [START_REF] Mermut | Baseline studies of The Clay Minerals Society Source Clays: Layer-charge determination and characteristics of those minerals containing 2 : 1 layers[END_REF] but lower than the value reported by [START_REF] Jaynes | Charge reduction, octahedral charge, and lithium retention in heated, Li-satured smectites[END_REF]. The mean layer charge shows an increase from 0.36 eq/(Si,Al) 4 O 10 after 5 days to a value of 0.42 eq/(Si,Al) 4 O 10 at the end of the experiment (Fig. 4a). The initial SAz-1 shows a heterogeneous charge distribution containing ~ 4 % of layers with a charge > 0.4 eq/(Si,Al) 4 O 10 (Fig. 4b). With increasing run time the layer charge distribution is shifted towards higher charge densities. After 150 days more than ~ 37 % of the smectite layers have a charge higher than 0.4 eq/(Si,Al) 4 O 10 (Fig. 4b).

XPS/SEM:

The Al/Si and the Mg/Si atomic ratios of the Na-saturated solids were measured to follow changes in smectite chemistry as a function of reaction progress. Both ratios increase as a function of reaction time indicating the increase of Al and Mg in the smectite structure. The rise of the Al/Si atomic ratio as function of a reaction time from 0.39 to 0.48 after 150 days was observed (Fig. 5). The Mg/Si ratio increased from 0.17 to 0.3 after 150 days (Fig. 5). The initial values of 0.39 and 0.17 measured on the initial smectite are in perfect agreement with the values measured by [START_REF] Metz | Stoichiometry of smectite dissolution[END_REF] for SAz-1 montmorillonite. SEM observations were carried out to verify the formation of precipitated aluminium or magnesium bearing phases which are below the detection limit of XRD. Some rare small brucite crystals but no Al-bearing phases were found after 80 days of reaction time with SEM. Although the information depth of the XPS signal is 6.3 nm, it cannot be concluded a priori that using the XPS only a surface layer of around five smectite layers is analyzed. Homogenously distributed contaminants may be detected by the XPS signal as well, because the analyzed sample area is about 1 mm 2 . However, accessory minerals which are within a smectite aggregate or covered by a smectite coating would not be included in the XPS analysis. Some SrCl 2 impurities were found by SEM and during XPS analysis.

STA measurements:

For the STA measurements we used sodium saturated samples. Based on the dehydroxylation temperature of 643°C the SAz-1 is a cis-vacant (cv) smectite [START_REF] Drits | An improved model for structural transformations of heat-treated aluminous dioctahedral 2:1 layer silicates[END_REF][START_REF] Emmerich | The cis-and trans-vacant variety of a montmorillonite: an attempt to create a model smectite[END_REF]. With increasing reaction time the dehydroxylation temperature decreases continuously to 605°C after 150 days (Fig. 6). In addition a second (small) dehydroxylation peak occurred at about 790°C after 80 days of reaction time which may be attributed to a trioctahedral smectite.

The dehydration temperature of the initial sodium saturated smectite was measured to be 147°C. For the first 80 days the dehydration temperature varies within the error of the measurement. After 150 days the dehydration reaction splits and a second maximum at 127°C was observed in DSC curves (Fig. 6). TG measurements demonstrate that the amount of interlayer water increased from about 12 % of the untreated Na-saturated SAz-1 to about 17 % after 150 days of reaction.

Infra-Red: IR spectroscopy revealed unchanged smectite spectra for all the experiments. No loss of intensity of vibration bands as a function of reaction time was observed. A weak band at 792 cm -1 indicating traces of poorly crystalline silica disappeared with increasing reaction time (Madejova and Komadel, 2001). No band at ~ 3680 cm -1 was detected in the reaction products which would indicate the presence of a trioctahedral smectite.

AFM: An AFM image of Na-saturated smectite particles after 150 days reaction time is shown in Fig. 7. The AFM analysis of the initial and the reacted particles indicate no major morphological changes throughout the experiment. Therefore, the following description includes all samples. The particles consist of irregular and fine grained plates. Mostly the edges are rounded, but straight edges appear as well. The top view shape of the particles resembles an ellipse. The height and the surface area of each of the selected particles were measured. In all 400 heights and 350 surface area measurements were performed. The distribution of the heights measurements for the initial smectite and the Reactor 6: 150 days sample show no change in particle morphology, particle height and surface area. The number of layers in the smectite particles was calculated assuming that the thickness of a single TOT layer with the hydrated sodium layer was about 12.5 Å [START_REF] Jasmund | Tonminerale und Tone. Struktur, Eigenschaften[END_REF]. Since diluted smectite suspensions were used and the samples were dispersed during the preparation for the AFM measurements the majority of the particles consists of 1-5 TOT layers (Fig. 7).

Oxygen Isotope Analyses

Discussion

The major objective of the present study is to demonstrate changes in mineral structure of a smectite at low temperatures under conditions of high alkalinity. Other studies have been performed on similar material but not a low temperature. Several different methods were used to obtain information on mineral evolution during the experiments.

Direct observation of crystal morphology

AFM observations of crystallite size, shape and thickness supports the hypothesis that the process is dominated by SST [START_REF] Baronnet | Silicate microstructures at the sub-atomic scale[END_REF]. The materials appear to retain their outer morphology but change their internal crystallographic composition.

Measurements of structural change

XRD

In contrast to [START_REF] Whitney | Dioctahedral smectite reactions at elevated temperatures: Effects of Kavailability, Na/K ratio and ionic strength[END_REF][START_REF] Whitney | Experimental investigation of the smectite to illite reaction: Dual reaction mechanisms and oxygen-isotope systematics[END_REF]Bauer and[START_REF] Bauer | Smectite transformation in high molar KOH solutions[END_REF] the initial and reacted smectites were measured in the Sr-saturated state at a fixed relative humidity. Under these conditions misinterpretation of dehydrated layers for illite layers can be neglected. XRD and infrared spectra gave little evidence for a change in smectite mineralogy in our experiments.

The number of expandable layers with one water layer increases, whereas the portion of the 0 water layers is about constant with reaction time. This suggests an increase in overall layer charge on the smectite. The peak width of the (001) reflections increases with reaction time suggesting a decrease in average diffracting domain. The peak position of the (060) band increases suggesting a change in internal mineral composition and an increase in charge on the structures. The layer charge using XRD and alkyl-ammonium treatments was found to increase from 0.32 to 0.42 and a clear shift in the layer charge distribution towards high charge layers was observed. Small amounts of a trioctahedral mineral were apparent (from (060) peak measurements) which indicates an increase of magnesium in the smectite.

MS measurements

The dehydroxylation temperature decreased from 643° to 605°C. This shift confirms the increased layer charge of the smectite. A lower dehydroxylation temperature for smectites with higher layer charge has been observed by Schultz (1969) and [START_REF] Michot | Hydration and swelling of synthetic Na-saponites: Influence of layer charge[END_REF]. In smectites, the octahedral sheet tends to be cis-vacant while the illites have mainly trans-vacant disposition [START_REF] Drits | Semiquantitative determination of trans-vacant and cis-vacant 2 : 1 layers in illites and illite-smectites by thermal analysis and X-ray diffraction[END_REF]. It is evident that the reacted smectite has a different disposition of the hydroxyl groups but that the cis-vacant character of the octahedral sheet remained nearly unchanged throughout the experiment. A second but small dehydroxylation appeared at 790°C. A similar observation was made by [START_REF] Emmerich | The cis-and trans-vacant variety of a montmorillonite: an attempt to create a model smectite[END_REF] for a de-and rehydroxylated Na-saturated cv smectite. XRD measurements of the 060 reflex indicate the appearance of trioctahedral Mg-smectite layers which is consistent with the XPS data showing an increase of Mg in the smectite structure (Fig. 5) and XRD data for the (060) peak indicating a small amount of a trioctahedral phase.

MS curves for the Na-saturated samples SAz-1 and the resaturated Reactor 6: 150 days sample are shown in Fig. 6. The amount of adsorbed water depends on the nature of the interlayer cations, the pre-treatment of the sample and the silica lattice [START_REF] Grim | Clay Mineralogy[END_REF]. Two types of the water molecules can be distinguished in the interlayer space of smectites: filling water which belongs to outer hydration shells of the cations or is bound to the surface of the smectite layers; and water belonging to the first hydration shell of the exchangeable cations.

Both smectites in Fig. 6 have the same interlayer cation and the same pre-treatment so that the differences can be attributed to changes in the silica lattice and as a consequence in the number of interlayer cations. It is evident that with increasing layer charge more cations are needed to compensate the charge deficit. This was confirmed by TG measurements which indicate an increase of the interlayer water from 12 to 17 % after 150 days of reaction time corresponding to increase of about 8 to 10 water molecules per cation Na + . These numbers are within the range of the values for Na-saturated saponites dominated by 2W layers [START_REF] Michot | Hydration and swelling of synthetic Na-saponites: Influence of layer charge[END_REF]. This is an important observation because it means that the hydration sphere of the Na + cations is already sufficient to induce the major presence of 2W layers from the beginning. This means that additional H 2 O molecules are not strongly bound to interlayer cations. In turn they should be easier to remove. In our experiments the occurrence of a second dehydration reaction centred at 127°C after 150 days was observed (Fig. 6). This is in contradiction to what someone would expect when the amount of absorbed water increases but it is consistent with the water vapor adsorption-desorption isotherms reported by [START_REF] Michot | Hydration and swelling of synthetic Na-saponites: Influence of layer charge[END_REF] that show a transition towards higher hydration states at lower P/P 0 values. As we know meanwhile hydrated smectites consist of %2W, %1W, %0W layers that superimpose in one peak which is often not resolved in DSC or MS curves. In addition water adsorbed on clay surfaces contributes to observed dehydration.

Measurements of chemical change

Oxygen constitutes the structural skeleton of the reacting materials, the degree of isotopic resetting can be used to indicate the extent of rearrangement of the mineral structure at high temperatures [START_REF] Whitney | Experimental investigation of the smectite to illite reaction: Dual reaction mechanisms and oxygen-isotope systematics[END_REF]. The initial rapid shift in δ 18 O values of 2.2‰ to more negative values for the reacted smectites is in agreement with an exchange of structural surface oxygen with water/KOH oxygen. This is due to the negative δ 18 O value of the exchange fluid (-7.4‰). According to [START_REF] Kohn | Oxygen isotope geochemistry of the amphiboles: Isotope effects of cation substitutions in minerals[END_REF][START_REF] Kohn | Oxygen isotope geochemistry of the amphiboles: Isotope effects of cation substitutions in minerals[END_REF][START_REF] Sheppard | Stable isotope geochemistry of clay minerals[END_REF][START_REF] Sheppard | Stable isotope geochemistry of clay minerals[END_REF] for temperatures of 80°C a fractionation of around +19‰ for smectite -fluid exchange applies, which would yield an end value of +12‰ for a smectite in equilibrium with the fluid.

Since at a solid/liquid ratio of 1/80 the δ 18 O value of the liquid can be regarded as constant throughout experiment run time. A rapid initial shift of the oxygen isotope ratio followed by a decrease in reaction rate with time in the reacted smectites can be observed. The isotope data can be taken to indicate a rapid equilibration of the solids with the aqueous solution. This would indicate no significant recrystallization of the smectite material.

XPS data showed that the increase of aluminium in the structure (increase of the Al/Si ratio from 0.39 to 0.48) is responsible for the increase in the layer charge. This is consistent with the solution chemistry which indicates an increase in silica in solution and magnesium relative to aluminium.

Chemical measurements of dissolved species and oxygen isotopes of solids indicates that little dissolution occurred but significant change in the overall composition of the solids is seen through changes in solution composition.

Comparison with previous experiments

Our results are consistent with previous studies which were conducted at neutral pH and higher temperatures (100° and 250°C) with various sources of potassium. Howard and Roy 1985 found the formation of high-charged layers when montmorillonite reacted at 150° and 250°C. [START_REF] Eberl | Hydrothermal reactivity of smectite[END_REF][START_REF] Whitney | Dioctahedral smectite reactions at elevated temperatures: Effects of Kavailability, Na/K ratio and ionic strength[END_REF] found that the formation of high-charge layers is promoted by an increase of potassium in the system. [START_REF] Yamada | Isothermal treatments of regularly interstratified montmorillomite-beidellite at hydrothermal conditions[END_REF] observed the formation of a beidellite plus saponite assemblage after the experimental alteration of synthetic beidellite-montmorillonite mixed-layers at high temperatures. [START_REF] Beaufort | An experimental alteration of montmorillonite to a di plus trioctahedral smectite assemblage at 100 and 200 degrees C[END_REF] observed during montmorillonite alteration at 100° and 200°C the formation of trioctahedral Mg-smectite and high-charge beidellitic layers. The beidellite + saponite assemblage is also observed between 100°-200°C in active geothermal fields where montmorillonite seems to be no longer stable at temperatures > 100°C [START_REF] Beaufort | Les interstratifies I-S et C-S dans les champs geothermiques actifs: sont-ils comparables ceux des series diagenetiques? Bulls[END_REF].

Process of mineral change

Given the observations of crystal structure (diffracting domain size, a-b cell size dimensions, direct AFM observations of crystal dimensions) it appears that there is a change with reaction time in the smectite structure. This seems to be due to the substitution of Al and Mg for Si in the smectite which creates smaller coherent diffraction domains in the structure without necessitating a dissolution -re-crystallization of the silicate mineral.

The observations cited above indicate that there is a significant, gradual change in the mineral cation composition of the 2:1 layers of smectite in the presence of KOH solution. The smectite minerals become higher charged, and have smaller coherent diffracting domains.

This involves significant chemical changes within the mineral crystallites. Aluminium and magnesium substitutes for silicon to a significant degree. The question is then: How does this happen? Essentially two pathways can be proposed according to classical theories: dissolution/crystallization or solid state transition.

If dissolution/precipitation is not the mechanism of mineral chemical change, one must accept a mechanism of diffusive modification within the existing crystallites. This change occurs at low temperatures, 80°C. Such a conclusion is surprising given the assumed very low rate of ionic diffusion in silicate crystallites at low temperatures [START_REF] Kyser | Retrograde exchange of hydrogen isotopes between hydrous minerals and water at low temperatures[END_REF].

However, the experimental evidence shown here strongly rejects the hypothesis of dissolution and crystallization, the other alternative to mineral change in aqueous solutions.

The conclusion of diffusion mechanism of mineral change at low temperatures is surprising.

Yet if one considers the great complexity and range of solid solution common in low temperature clay minerals, one can suspect that important mineral compositional change occurs at these low temperatures. Normally one attributes high solid solutions in solids to high thermal energy but it is known that chemical variability can also increase at lower temperatures in systems such those where illite is stable [START_REF] Velde | Clays Minerals: A physico-chemical explanation of their occurence[END_REF].

We propose that the phenomenon of chemical change at low temperatures is one of a response to new chemical conditions inducing a larger range in chemical composition which, in the case studied, increases the alumina and magnesium content of the smectite minerals observed.

In this particular instance the change is accomplished by internal diffusion within the smectite crystallites. The relative contribution (S cont.) of the two structures S1 and S2 is given in % as well as the proportion of bi-hydrated, mono-hydrated and dehydrated layers (2W, 1W and 0W layers, respectively) within these elementary contributions. The total proportion of 2W, 1W and 0W layers in crystal is given in % as Tot. 2W, Tot. 1W and Tot. 0W, respectively. Layer thickness (LT) of the different layer type is given in Å. N is the mean number of layers in the coherent scattering domain.

Tables

The structure model of initial sample is given from the study of [START_REF] Ferrage | Etude expérimentale de l'hydratation des smectites par simulation des raies 00L de diffraction des rayons X. Implications pour l'étude d'une perturbation thermique sur la minéralogie de l'argilite du site Meuse-Haute Marne[END_REF] Table 2. Optimum structural parameters used for the simulation of XRD profiles.

:

  The oxygen-isotope composition of the Na-saturated clays was monitored to provide another measure of reaction progress and is shown for the initial and the reacted smectites in Fig. 8. The δ 18 O values of the clays changed rapidly from δ 18 O = 17.4‰ for the initial SAz-1 to a value of δ 18 O = 15.2 ‰ for the reacted smectite after 150 days. A similar trend was observed by Whitney and Northrop, 1988 at higher temperatures (250°-450°C) for the SWy-1 smectite in distilled water. The composition of the water in the experiment was -7.5 ‰. The isotope values change with time, decreasing rapidly in the beginning before attaining a value near 15.2‰, which shows only partial exchange with the aqueous solution oxygen. The initial exchange is most likely one of surface interaction. If dissolution -recrystallization occurred one would expect greater exchange between solution and silicate material.

Figures

  Figures

Figure 1

 1 Figure 1Comparison between experimental (crosses) and calculated (solid line) XRD patterns for the reacted Sr 2+ saturated SAz-1 recorded at 40% relative humidity. The grey bars outline the modifications in peak width between 5 days and 150 days of reaction time.

  Figure 2Evolution of the 060 bands as a function of reaction time.

Figure 5

 5 Figure 5

  Mg analyses was done using synthetic SiO 2 , Al 2 O 3 and MgAlO 4 samples. The SiO 2 (30 nm SiO 2 layer on silicon wafer), Al 2 O 3 (sapphire single crystal, TBC-Kelpin Co., Neuhausen) and

	H 2 O/O 20 (OH) 4 in 1W and 2W layers, respectively. Other structural parameters and atomic
	positions in the different layer types are described in Ferrage et al. (2005).
	Briefly, the fitting strategy developed by Ferrage et al. (2005) consists in using a main
	structure, periodic when possible, to reproduce as much as possible of the experimental XRD
	pattern. If necessary, additional contributions to the diffracted intensity are introduced to
	account for the misfit between calculated and experimental patterns. The presence of two
	MLSs does not imply that two populations of particles are physically present in the sample.
	As a result, layers with the same hydration state present in the different MLSs contributing to
	the diffracted intensity are assumed to have identical properties (chemical composition, layer
	thickness, z-coordinates of atoms). The relative proportions of the different MLSs and that of
	the different layer types in these MLSs were also considered as variable parameters during the
	fitting procedure.
	. In addition, smectite hydration heterogeneity was monitored using the
	methodology developed by Ferrage, 2005a. This heterogeneity results from the coexistence in Layer Charge: A detailed description of the experimental protocol and the layer charge
	the same crystal of the different layer types usually considered in smectite structure in non-determination can be found in Mermut and Lagaly, 2001 and Hofmann et al., 2004. All
	saturated water vapour conditions: dehydrated layers (0W layers, layer thickness at 9.6-measurements were done in triplicate.
	10.1 Å), mono-hydrated layers with one plane of H 2 O molecules in the interlayer (1W layers
	at 12.3-12.7 Å), and bi-hydrated layers with two planes of H 2 O molecules in the interlayer X-ray Photoelectron Spectroscopy (XPS): XPS analyses were performed using a Physical
	(2W layers at 15.1-15.8 Å). Such structural heterogeneity can be seen as the result of both Electronics Inc. PHI 5600ci spectrometer equipped with a hemispherical capacitor analyzer
	contrasted layer charges from one layer to the other or within a single layer and intrinsic
	thermodynamic effects.
	For the modelling of all experimental XRD pattern, the preferred orientation of the sample
	(σ*) was set to 7°. Because the water content of hydrated layers cannot be determined
	accurately from the sole 001 reflection, the number of H 2 O molecules was set to 3 and 6

and multichannel detector. Si 2s, Al 2p and Mg 2s elemental lines were recorded by Al K a Xray excitation (1486.6 eV) using a source power of 200 W, a source -sample distance of 10 mm and a source -analyzer angle of 54°. The angle between sample surface normal and analyzer was 45°. Atomic concentrations were calculated from peak areas of elemental lines after Shirley background subtraction, by use of sensitivity factors and the analyzer transmission function determined at same acquisition conditions. Calibration of the Si, Al and

Acknowledgement:

We thank Klaus Spieler and Eva Soballa (INE) for their technical assistance during the SEM measurements and Dieter Schild (INE) for performing the XPS analysis. We would like to thank Heike Pieper and Dirk Bosbach for their help performing the AFM measurements.