Waseem Al Hadad

Vincent Schick

Denis Maillet
email: denis.maillet@univ-lorraine.fr

Fouling detection in a shell and tube heat exchanger using variation of its thermal impulse responses : methodological approach and numerical verification

Keywords: Heat exchanger, fouling detection, impulse response, transfer function, deconvolution, conjugate heat transfer. Nomenclature a thermal diffusivity, m 2

A new method for fouling detection in a heat exchanger is proposed. This method is based on the sensitivity of its impulse response (time version of its transfer function) with fouling. The concept is validated here numerically for a shell and tube heat exchanger in steady state regime, perturbed thermally by a transient surface source. This method is presented here, on a theoretical ground. It can be implemented online, whatever the geometric and structural characteristics of the setup. It only requires the use of regular temperature sensors and actuators. The only condition concern the equations of the model of the heat exchanger which have to be linear with time invariant coefficients.

x, y, z space coordinates, m

Introduction

During the operation of a heat exchanger, deposits may develop over the wall that separates the hot and cold fluids. Over time, these materials, which are either layers sticked to the walls or porous materials in the volume of the flow, appear and grow along the length of wall or/and in specific areas. They can thus block or change the flow direction. Furthermore, they are less conductive than the separating. Therefore, they will significantly reduce the heat transfer efficiency and performance in the course of time.

In the literature, this undesirable fouling layer degrades the operation of the heat exchanger and one of the major difficulties is the early quantification (prediction) or detection of this undesirable layer. Fouling quantification is always more difficult than fouling detection. This difficulty comes from the random behavior of this phenomenon.

There are various methods to slow down the presence and accumulation of fouling (filtering, chemical treatment, ..) [START_REF] Bott | Fouling of heat exchangers[END_REF][START_REF] Diaz-Bejarano | Organic and inorganic fouling in heat exchangers: Industrial case study analysis of fouling rate[END_REF]. As soon as the fouling takes place, the complete cleaning of the installation becomes essential. The systematic periodic cleaning to remove fouling is not the best solution. If the period is short, the costs of production shutdown, that is production losses during cleaning and maintenance tasks (cleaning, disassembly,...), can increase. On the contrary, if the this period is long, fouling can cause an increase of the energy consumed by the process (the heat exchanger efficiency decreases with fouling), an increase in the consumption of cooling fluid and an impact on the environment. For all these reasons, the manufacturers of heat exchangers and their customers seek to move thus from systematic maintenance to predictive maintenance.

Apart from the methods based on models of the black box type, such as neural networks [START_REF] Guo | Optimized neural network-based fault diagnosis strategy for VRF system in heating mode using data mining[END_REF], popular methods of fouling detection are based on the decrease in the Overall Heat Transfer Coefficient, noted U. In some configurations, this method is not precise enough. This can be explained by the fact that the models used for explaining U are based on a strong assumption (locally 1D transfer which consists in neglecting the axial heat transfer in the solid walls) and by the compensation that can take place between the resistances to conduction and advection : the fouling layer tends to increase the resistance to conduction (due to the increase of the wall thickness) and to decrease the resistance to convective heat transfer (with fouling, the hydraulic diameter decreases, which increases the internal velocity if the flow rate is imposed and that will lead to an increase of convective heat transfer coefficient). Therefore, the overall exchange coefficient does not always have a good sensitivity to fouling [START_REF] Hadad | French : Thermique des mini-canaux : comportement instationnaire et approche convolutive)[END_REF]. Another classical method of detection, easy to implement, is based on the increase in pressure loss when the flow rate is imposed but the effectiveness of this method depends on the type of exchanger [START_REF] Gudmundsson | Comparison of fouling detection methods using experimental data[END_REF].

Improved results can be reached if one uses transient regimes, such as a steep variation of one of the flowrates, which allows a better followed up of the Overal heat transfer coefficient, using dynamic models and state estimation through Kalman filterng, see [START_REF] Gudmundsson | Online fouling detection of domestic hot water heat exchangers[END_REF][START_REF] Palmer | Analysis of transient data in test designs for active fault detection and identification[END_REF] or variation of the fluid inlet temperatures [START_REF] Pham | Fouling detection in a parallel flow heat exchanger via a Roesser model identification procedure[END_REF].

Another possibility is to work in modulated regime, using perturbations of the flowrates or inlet temperatures [START_REF] Lalot | Analysis of the unsteady behavior of an electrical heater using the lock-in technique: Comparison of theory and experiments[END_REF][START_REF] Lalot | The lock-in technique applied to heat exchangers: A semi-analytical approach and its application to fouling detection[END_REF], which allows the acquisition of signals of good signal to noise ratios through lock-in techniques.

A common approach for fouling detection consists in comparing experimental temperatures with those predicted by a reduced model developed earlier [START_REF] Jonsson | Use of extended Kalman filtering in detecting fouling in heat exchangers[END_REF], [START_REF] Delrot | Fouling detection in a heat exchanger by observer of Takagi-Sugeno type for systems with unknown polynomial inputs[END_REF]. Recently Palmer et al. [START_REF] Palmer | Optimal design of tests for heat exchanger fouling identification[END_REF][START_REF] Palmer | Analysis of transient data in test designs for active fault detection and identification[END_REF] attempted to quantify fouling in a cross-flow plate fin heat exchanger. The idea was based on the minimization of the differences between the measured temperatures and their predictions by a model that takes fouling into account. We point out here the common simplifying assumptions in these works : the fluid flow was considered one-dimensional along each fluid direction, the axial heat conduction along the fluid was considered negligible and the convective heat transfer coefficient was considered uniform along each fluid direction and the heat transfer between the two fluids was locally 1D. The conduction in the walls has been neglected in [START_REF] Jonsson | Use of extended Kalman filtering in detecting fouling in heat exchangers[END_REF] and [START_REF] Delrot | Fouling detection in a heat exchanger by observer of Takagi-Sugeno type for systems with unknown polynomial inputs[END_REF] and considered 2D in [START_REF] Palmer | Optimal design of tests for heat exchanger fouling identification[END_REF] (with a zero temperature gradient in the direction of the thickness).

It should be emphasized that in the case where the Péclet number Pe is small with non negligible thickness of the walls and short channels, the axial conduction in these walls plays a very important role [START_REF] Herwig | Flow and heat transfer in micro systems : Is everything different or just smaller ?[END_REF][START_REF] Herwig | Critical view on "new results in microfluid mechanics : an example[END_REF][START_REF] Guo | Size effect on microscale single-phase flow and heat transfer[END_REF]. So, under these circumstances, these methods do not provide a good detection of fouling.

There are also sophisticated methods that have been developed to detect fouling, such as silicon sensors [START_REF] Stenberg | A silicon sensor for measurement of liquid flow and thickness of fouling biofilms[END_REF], ultrasonic, acoustic and optical techniques [START_REF] Withers | Ultrasonic, acoustic and optical techniques for the noninvasive detection of fouling in food processing equipment[END_REF]. These methods which require a specific instrumentation, are not always adapted to the equipment and are therefore limited to a laboratory use.

So, the challenge is to be able to detect fouling in real time, without stopping the process and implementable at both the laboratory and industrial scales.

In this work, we propose a detection technique based also on comparison of experimental measurements and outputs of a reference model, that is identifying the impulse response of the system. The use of this concept of impulse response is based on the only assumption that the system of equations that model the response of the system to a perturbation is linear, with time invariant coefficients (LTI). This means that this type of method could allow getting rid of the other a priori arbitrary asumptions required by the models used in the detection methods analyzed above, which could add some robustness to it.

So in the detection technique proposed here, we compare the system signature obtained experimentally during operation with the reference signature. This reference model must be identified previously from a calibration experiment when the exchanger is clean. This signature or model has a convolutive nature, which is valid for a LTI system. The system signature corresponds here to its impulse responses or its transfer functions (the transfer function is the Laplace transform of the impulse response).

These functions, which can be identified whatever the geometric complexity (conjugated heat transfer, turbulent, flow direction change, roughness, development of the boundary layer,..), depend only on the geometry, of the velocity field of the flowing fluids and of the correponding thermo-physical properties. These features can therefore be used to detect the temporal variation of a structural parameter of the system during operation. In our previous works, we have shown the interest of transfer functions for LTI systems where heat diffusion and advection occur in a conjugated way, both on a theoretical [START_REF] Hadad | Modeling unsteady diffusive and advective heat transfer for linear dynamical systems : A transfer function approach[END_REF] and on an experimental [START_REF] Hadad | Experimental transfer functions identification : Thermal impedance and transmittance in a channel heated by an upstream unsteady volumetric heat source[END_REF][START_REF] Hadad | Inverse conduction and advection in a flat channel with transient external thermal excitation and observation[END_REF] basis. These are either impedances relating a local or averaged temperature to the cause of its variation, that is an unsteady heat source, or a "transmittance" linking a local or averaged temperature to another one whose variation occurs before. In that last case, the prior variating local temperature will be called a "pseudo-source" here, since it is also the consequence of the common heat source.

Our paper is organized as follows : in section 2, we show how the transfer function can be identified whatever the geometry and how it can be used to detect the possible presence of fouling. In section 3, we present the studied system (shell and tube heat exchanger) used in this work to validate the thermal fouling detector. In section 4, we apply this methodology to the model of the above system, using Finite Element simulations [START_REF]COMSOL Multiphysics version 5.3[END_REF].

Impulse response identification and principle of a thermal fouling detector

This section presents the general concept of identification of the impulse response of a linear and time invariant system (LTI), using the knowledge of both input (the perturbation) and output (the corresponding response).

Once an impulse response has been identified for such a system (a calibration), its re-identification at a later time, with comparison with its prior values, gives information about a possible change of one of its structural parameter (a non destructive test). This principle will be applied to a shell and tube heat exchanger in the next chapters.

In one of our works just cited above [START_REF] Hadad | Modeling unsteady diffusive and advective heat transfer for linear dynamical systems : A transfer function approach[END_REF], we focused on a physical system modeled by a LTI mathematical system composed of a partial differential equation, such as the heat equation and of its associated boundary and interface conditions. We showed that, for such a system, the output (response or consequence) y at a given point, or its average in any subsurface or subdomain, is a simple product in Laplace domain, or a convolution product in time domain, between the input (excitation or cause) u and a corresponding function

H : ( ) ( ) ( ) (Laplace domain) = y p H p u p (1a) ( ) (t) ( ) (time domain) = y t H * u t (1b)
where the upper bar ( _ ) denotes the Laplace transform and p the Laplace parameter : 0 ( ) ( ) exp(-)d for or

ψ ψ ψ ∞ = = ∫ p t pt t y, H u (1c) Function (t)
H is the impulse response, while its Laplace transform ( ) H p is the corresponding transfer function. The star (*) denotes the convolution product.

The convolution product (1b) can be defined in two different ways, using its commutative property [START_REF] Chata | Estimation of an aerosol source in forced ventilation through prior identification of a convolutive model[END_REF] : ). This can be put under a column vector/matrix form:

0 ( ) ( ) ( ) d = - ∫ t y t H t t' u t' t' (2a) or 0 ( ) ( -) ( ) d t y t u t t' H t t' = ∫ (2b)
= ( ) y H u M (3a)
We use the commutative property of the convolution product here :

= ( ) y u H M (3b) with 1 1 2 1 2 3 2 1 3 1 2 1 0 0 0 0 0 ( ) = with t t t t N N N N t ψ ψ ψ ψ ψ ψ ψ ψ ψ ψ ψ ψ ψ ψ - -                 ∆ =                 ψ ψ L O M O M M M O M M K M (4) and
( )

1 1 1 1 ( ) d ( ) + ( ) for or 2 k k t k k k t t t t t u H t ψ ψ ψ ψ ψ - - = ≈ = ∆ ∫ (5)
with 0 0 t = , 0 ( ) = 0 t ψ and 1 to t k N =
. Here (.) M is a function whose value is a square matrix of size x t t N N (a lower triangular Toeplitz matrix) that depends on a column vector, either H or u .

Let us note that the numerical quadrature, in equations (3a) and (3b), of the convolution products (2a) and (2b) respectively, is only valid for a time step small enough with respect to the chararactistic times of both input ( ) u t and impulse response . We also assume that these two functions are continuous and bounded over the 0 final , t     interval. If there is a discontinuity at time 0 0 t = (case met if ( ) t ψ is a step function) , the first coefficient of vector ψ gets the following value :

1 1 ( ) t ψ ψ = .
So, if both input and output are known, the impulse response H can be identified by inverting the Toeplitz matrix of the source, that is ( ) u M , see equation (3b).

Once the impulse response has been identified for a given system, it can be considered as its identity card : the impulse response is an intrinsic quantity that depends only on the structural characteristics of the system : dimensions (wall thickness, length, diameter, etc.), thermophysical properties and internal velocity field.

During the system operation, if the impulse response has changed, this can be explained by a change in the structural parameters of the system.

In the case of a heat exchanger, such a change may stem from a change of one of the fluid mass flow rates or from the apparition of fouling. On this basis, a thermal detection method can be designed for fouling. Its principle can be summarized in the following steps :

1. The type of thermal perturbation ( ) u t and the location of its application have to be chosen, together with the place the output ( ) y t is obseved. This perturbation should be small enough in order for the response of the system to stay in the linear domain.

The impulse response ( )

clean H t is identified when the system is clean or new, for given mass flow rates of the fluids. This principle of fouling detection will be verified numerically in the next sections for a tube and hell heat exchanger.

Remark : For some specific temporal forms of the source ( ) u t the transfer function and the corresponding impulse response, can be calculated without using the solution of the vector/matrix linear system (3b). This alternate method, which is semi-analytic, is possible only for a temporal source form corresponding to available explicit inverse Laplace transforms of both source and response.

Here, we consider the excitations of the following forms: a Dirac pulse, a step and a ramp. For these special excitation forms ( ) We also assume that the heat exchanger initially operates in steady state and at a given moment, the transfer is perturbed by a surface heat source, ( ) q t which is uniformly distributed over the outer lateral surface of the inlet hot fluid pipe over a length 3 l , see Figure 1. In this application, ( ) q t corresponds to the input denoted by ( ) u t in the general presentation in section 2.

u
The two volume flow rates are assumed to remain constant in time during the thermal perturbation and both fluid and solid walls have thermophysical properties that are supposed to be independent of temperature : the equations of motion can be solved independently of the heat equation. The flow in the tube is assumed to be laminar (Re < 2000) and is described by the classical conservative equations of mass and momentum.

Turbulence in the shell is taken into account through the k ε

model [START_REF] Jones | The prediction of laminarization with a two equation model of turbulence[END_REF] :

. = 0 ∇ v ( . ) = ( + ) ( ( ) )

T T P . ρ µ µ   ∇ -∇ + ∇ ∇ ∇ +   F v v v + v (7)
( . ) = . ( )

T k k k k P µ ρ µ ρ ε σ   ∇ ∇ + ∇ + -     v (8) 2 1 2
( . ) = . ( )

T k C P C k k ε ε ε µ ε ε ρ ε µ ε ρ σ   ∇ ∇ + ∇ + -     v (9) with 2 T C k / µ µ ρ ε = , : ( ( ) ) T k T P µ   = ∇ ∇ ∇   v v + v , 1 1 44 C . ε = , 2 1 92 C . ε = , 0 09 C . µ = , 1 k σ = and 1 3 . ε σ =
. v is the velocity vector, the density, ρ the density, µ the dynamic viscosity, T µ the turbulent viscosity, F the external force per unit volume, k the turbulent kinetic energy, ε the rate of dissipation of the turbulent kinetic energy, while k

σ , ε σ , 1 C ε and 2 C ε are empirical constants.
A change of function is made in order to define the temperature variation θ , that is the response to the surface heat perturbation ( )

q t : ( ) = ( ) ( 0) x, y, z, t T x, y, z, t T x, y, z, t θ - = (10)
The transient heat equation in the walls (tube, shell and baffles walls), in the fluid flow (hot and cold), the initial conditions and the boundary conditions (for heat and motion equations), are :

• Heat equation in the walls and in the fluid flows :

2 1 1 . i i i i a a t θ θ θ ∂ ∇ - ∇ = ∂ v with or i s f = (11)
Here, subscript f designates the fluid domain while s corresponds to the solid part of the system, where = 0 v .

• Initial condition ( 0) 0 i x, y,z, t θ = = (12)
This condition means that the thermal perturbation ( ) q t is applied at time 0 t = on a system where a steady state had be reached for the initial temperature field ( 0) T x, y, z, t = .

• Boundary conditions -at the fluid inlets, temperature variations and volume flow rates are imposed: Inlet hot fluid temperature ( )

h f in t θ θ = (13a)
Inlet cold fluid temperature ( )

c f in t θ θ = (13b)
Inlet hot fluid volume flow rate

h v v Q Q = & & (13c)
Inlet cold fluid flow rate

c v v Q Q = & & (13d)
-at the fluid outlets, heat flux and pressure are imposed.

Diffusion is neglected with respect to advection, because the Péclet numbers Pe (based on the outlet diameters 1 D and 2 D ) are high in our application :

. 0 flux -= n q [START_REF] Herwig | Flow and heat transfer in micro systems : Is everything different or just smaller ?[END_REF] Pressure is also assumed to be uniform at the outlet of both hot and cold fluids, with a zero level here :

0 P = (15)
-the external fluid surfaces of the tubes and shell are supposed to be insulated : . 0 flux -= n q [START_REF] Guo | Size effect on microscale single-phase flow and heat transfer[END_REF] Equation ( 16) is valid everywhere, except over the cylindrical surface of length 3 l at the inlet of the hot fluid, see Figure 1, where the surface heat source ( ) q t is applied ;

. ( )

flux q t - = n q (17)
-at the solid/fluid interfaces, a continuity of temperatures and heat fluxes is assumed, with a zero local velocity.

In this work, the velocity field is assumed stationary and the thermophysical properties in the solid and in the fluid remain independent of temperature.

As a consequence, the considered system of equations (heat conduction and advection) ( 11) -( 17)) that governs the transient temperature field θ is linear and time invariant, since the coefficients of these equations (including the velocity field v ) do not depend on time.

So, the convolution model (Equation 1b, in its scalar form, where y is a temperature difference and u is the perturbation heat flux q, remains valid. In this case, the temporal transfer function (impulse response) H is called an impedance is denoted here Z , see [START_REF] Hadad | Modeling unsteady diffusive and advective heat transfer for linear dynamical systems : A transfer function approach[END_REF]. The parameterized form = Z H of the impedance can thus be identified by inverting vector equation (3b).

Results and discussion

The objective of this part is to validate numerically the principle of thermal fouling detector presented in section 2. The system represented in section 3 (tube and shell heat exchanger) has been modeled in 3D transient and simulated (after a mesh sensitivity study) by the commercial code COMSOL [START_REF]COMSOL Multiphysics version 5.3[END_REF].

Two numerical experiments were made with two different temporal forms of the surface heat source ( ) q t (thermal perturbation). In the first experiment, a constant surface heat source was applied from 0 to final t t t = = (step function) where final t is the experiment duration. In the second experiment, a surface heat source that varies linearly with time, 0 to final t t t = = (ramp function) was applied, such as ( ) = q t K t where K is the slope and t is the time.

In these two experiments, we assumed that the flowing fluid in the tube (hot fluid) as well as that in the shell (cold fluid) were water and the solid walls (tube, shell and baffles) were supposed to be made out of steel (1C).

The corresponding thermophysical properties, that were assumed to be independent of temperature, are given in Table 2.

Table 2

Thermophysical properties of each medium.

f λ (W.m -1 .K -1 ) f c ρ (kJ.m -3 .K -1 ) s λ (W.m -1 .K -1 ) s c ρ (kJ.m -3 .K -1 ) 0.63 4186 44 3666
The geometric parameters of these simulations, see Figure 1 are given in Table 3.

Table 3

Geometric parameters. Symbols The average velocity, m U calculated at the inlets, the volume flow rate, v Q & and the corresponding Reynolds number, Re and Péclet number, P e are given in Table 4 : here the same volume flow rates were chosen for the hot and cold fluids,

h c v v v Q Q Q = = & & & .
Since the inlet diameters were equal for the two fluids ( 13

D D = ,
see Figure 1 and Table 3), the average velocities were also equal

Table 4

Average velocity and corresponding dimensionless numbers.

1 (m.s ) m U - 3 1 (m .s ) v Q - & 1 m f Re U D / υ = P e Re Pr =
Hot and cold fluids 0.01 1.54 10 -6 140 980

Initially, in both numerical experiments, before thermal excitation (for 0 t ≤ ), the temperature fields in the solid and in the fluid were stationary, so ( ( 0) = 0 x, y, z, t θ = ), see equations ( 10) and ( 12), and both inlet temperatures remained constant during the thermal excitation. This means that the only source of heat was the surface heat source ( ) q t . The temperature responses at each point in the system correponding to these two excitations forms ( ) q t were calculated at each instant i t for 1 to The impulse response H that connects the excitation (the thermal or power input) to its response in temperature, is noted Z . We consider here that the source is q and its responses, the average temperatures The evolution of the heat source, for each experiment (step and ramp) and of its responses are plotted in Figure 3 (step) and in Figure 4 (ramp). Starting from these synthetic profiles, the impulse responses, h Z and c Z , have been identified using the convolutive model (Equation 3b) or by the semi-analytic expressions given in Table 1. These identified impedances are plotted in Figure 5, for step and ramp excitation.

t i N = with t final N t / t = ∆ , for a time step 1 s t ∆ = .
The initial zero values for h Z and c Z , (for 9 s t <

) mean that temperature signals h out θ and c out θ need a finite amount of time to react (the thermal disturbance has not yet arrived). For large times, the impulse response returns asymptotically to zero at 250 s t  for c Z and earlier for h Z . This corresponds to the maximum duration of the response and can be considered as a local relaxation time of the system or a time necessary to reach a final steady state level in the outlet cross sections for a step excitation. In between these two extreme times, a unique maximum exists that gives information about the level of each response. It will be shown in a remark in subsection 4.3 further on that the larger the area below each impulse response is, the larger the corresponding "resistance" between the source and the outlet temperature of each fluid becomes.

Figure 5 shows that the modeled shell and tube heat exchanger behaves as a LTI system, since the identified impulse response, here an impedance, is independent of the excitation form (it is an intrinsic quantity and can be considered as a signature of the system). It also shows that the two methods of identification give the same results. 1) and * (identified from ramp and by Table 1). In this section, we will show the sensitivity of the impulse response to fouling. A fouling layer (limestone deposit) of uniform thickness c e is added over the hot face of the wall separating the two fluids in the system considered above, see Figure 6. Z . These figures show that the sensitivity of the impulse responses to fouling are high (higher than the corresponding impedance for outlet fluid observation in the hot fluid case, see Figure 7) and that they can be used non-intrusively to detect fouling in the present case.

Impulse response identification from h

Interest of integrating the identified impedances with respect to time

As it was shown in [START_REF] Hadad | Modeling unsteady diffusive and advective heat transfer for linear dynamical systems : A transfer function approach[END_REF][START_REF] Hadad | Experimental transfer functions identification : Thermal impedance and transmittance in a channel heated by an upstream unsteady volumetric heat source[END_REF], an impulse response or impedance ( ) Z t can be written in steady state and this writing corresponds in fact to a thermal resistance, noted here ss Z , superscript ss designating a steady state value here. The latter can be calculated from the impulse response itself, see equation (18a) or from the asymptotic value of the source ss q and of the corresponding response ss θ , see equation (18b) :

0 ( ) d ss Z Z t t ∞ = ∫ (18a) ss ss ss Z / q θ = (18b)
By integrating each impulse response profile (see Figure 7, Figure 8, Figure 10 and Figure 11) for step or ramp excitation or by calculating its steady state value by the ratio of the asymptotic value of excitation ss q to the response ss θ (these last values are available for the step excitation only) one obtains the corresponding steady state impedance (it corresponds to a generalized thermal resistance), ss Z . The obtained values are given in Table 5.

Table 5

Identified steady state impedances (in K.m 2 .W -1 ) from the temporal distribution of the tranisent impedances (impulse response) or from the asymptotic values of the source ss q and the responses ss θ . Z for c e = 2 mm, because in this case, presence of the fouling insulating layer favours direct transfer from the source to the cold fluid at the expense of the hot fluid.

Internal measurements

Variation of steady state impedances, that are generalized thermal resistances, constitutes a convenient tool for fouling detection, even it they have to be interpreted with some hindsight. Their transient versions bring many more informations since they are functions and not scalars [START_REF] Hadad | Characterization of a heat exchanger by virtual temperature sensors based on identified transfer functions[END_REF].

Conclusions

We have proposed here a method for fouling detection in a heat exchanger or any thermal device whose model equations are linear and time invariant. The presented method is based on the variation of the impulse response that characterizes the system (thermal signature) with fouling. The idea here is to identify the model that describes the transfer from measurements on the system itself, that is via an experimental calibration (the identification phase). This identification, on a clean and then on an operating system, can replace fouling detection by a classical method based on estimation on the thermal parameters of a detailed model, that is the equations of the system, with boundary conditions and initial conditions with simplifying assumptions. This allows to make assumptions that are unavoidable in a classical modeling attempt. Indeed, the identified impulse response takes into account all the effects that are difficult to model implicitely: wall roughness, turbulence, contact resistance, defect in the materials that can have thermophysical properties different from their a priori known nominal values, etc.

In order to validate the concept of this method, it was applied here to output data of a numerical model of a tube and shell heat exchanger, where a thermal disturbance was made by a surface heat source located upstream the hot fluid inlet. The temperature responses were calculated at the outlets of the hot and cold fluids or on the external wall surfaces of the corresponding cross sections using a Computational Fluid Dynamics code, COMSOL. Three results were obtained :

i)

The impulse response, that depends on the definition of the way the thermal excitation is applied and on the location where the temperature response is measured, is able to reproduce the temperature outputs of the detailed model through its use in a specific reduced model, that is the convolution product with the input ; ii)

The estimated impulse response does not depend on the time shape of the excitation ; iii)

The impulse response is sensible to the presence of a fouling deposit.

So, this concept verification study shows that the impulse response is an intrinsic quantity and can be used to detect fouling in a heat exchanger. Let us note that the present technique has already been tested experimentally for a flat plate heat exchanger, see [START_REF] Hadad | Fouling detection in a heat exchanger using variation of its thermal transfer functions : an experimental approach[END_REF].

1 v

 1 coefficient, W.m -2 .K -1 m U mean velocity, m.s -velocity field, m.s -1

3 .Figure 1 .

 31 Figure 1. We assume that the circulating fluids in the tube (hot fluid) and in the shell (cold fluid) are incompressible.

Fig. 1 .

 1 Fig. 1. Shell and tube heat exchanger and correponding dimensions.

  the two fluids were the same, therefore the correponding dimensionless numbers were equal (

4. 1 .

 1 Impulse response identification from h out θ and c out θ without fouling(calibration)

  at the outlet of the hot and the cold fluids are h out θ and c out θ . The corresponding impulse responses are h Z and c Z respectively, see Figure 2.

Fig. 2 .

 2 Fig. 2. Shell and tube heat exchanger and correponding impedances h Z and c Z .

Fig. 3 .

 3 Fig. 3. Step excitation q and corresponding average hot and cold outlet temperatures ( h out θ and c out θ ).

Fig. 4 .

 4 Fig. 4. Ramp excitation q and corresponding average hot and cold outlet temperatures ( h out θ and c out θ ).

Fig. 5 .

 5 Fig. 5. Identified h Z and c Z from step and ramp excitations : continuous line (identified from step and equation (3b)), o (identified from ramp and equation (3b)), + (identified from step and by Table1) and * (identified from ramp and by Table1).

Fig. 6 .

 6 Fig. 6. Shell and tube heat exchanger with fouling.

Fig. 7 .

 7 Fig. 7. Identified h Z without ( c e = 0) and with fouling layer of uniform thickness c e = 1mm and c e = 2mm.

Fig. 8 .

 8 Fig. 8. Identified c Z without ( c e = 0) and with fouling layer of uniform thickness c e = 1mm and c e = 2mm.

4. 3 .Fig. 9 .

 39 Fig. 9. Shell and tube heat exchanger and correponding impedances w h Z and w c Z .

Fig. 10 .

 10 Fig. 10. Identified impedances w h Z without and with a fouling layer of uniform thickness c e = 1mm and c e = 2mm.

Fig. 11 .

 11 Fig. 11. Identified impedances w c Z without and with a fouling layer of uniform thickness c e = 1mm and c e = 2mm.

Table 1 Alternate

 1 t and from equation (1a), one can calculate semi-analytically the corresponding temporal transfer function ( ) H t using the corresponding responses, see Table1. These expressions of ( ) H t are valid only in cases where the response is not instantaneous, that is

	d ( 0 ) = 0 d y t t

+ =

. method for identifying the impulse response ( ) H t starting from responses ( ) y t to specific source forms ( )

u t ss u / p 2 K / p

Acknowledgment

The authors want to thank the FEDER European Fund as well as Région Grand Est who funded this work run within the EPHAISTOS project.


Eleanor Bakker

Bruno Lanson
email: bruno.lanson@univ-grenoble-alpes.fr

Nathaniel Findling

Michelle M Wander

Fabien Hubert

Mineralogical differences in a temperate cultivated soil arising from different agronomic processes and plant K-uptake

Keywords: Potassium (K), plant nutrition, fertilisation, clay mineralogy, X-ray diffraction, long-term field experiment, Morrow Plots

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Potassium (K) is an essential plant nutrient, playing roles in photosynthesis, water regulation and plant resistance to diseases [START_REF] Marschner | Mineral nutrition of higher plants[END_REF]Mengel et al., 2001). In soils, 98-99 % of the potassium is found in the crystal structure of K-bearing minerals such as Kfeldspars and micas (muscovite and biotite [START_REF] Brouder | Potassium Cycling[END_REF]Mengel et al., 2001;[START_REF] Sparks | Potassium dynamics in soils[END_REF]. In addition to climate, vegetation and factors like soil texture and pH, content and composition of soil solution, the nature and extent of weathering of K-bearing minerals present control the release of mineral K into solution through dissolution [START_REF] Manning | Testing the ability of plants to access potassium from framework silicate minerals[END_REF]Mengel et al., 2001;[START_REF] Sparks | Potassium dynamics in soils[END_REF][START_REF] Velde | Soils, plants and clay minerals[END_REF]. Apart from this mineral pool, K is found also as an exchangeable cation (1-2 %) and in the soil solution (0.1-0.2 % - [START_REF] Brouder | Potassium Cycling[END_REF]Mengel et al., 2001;[START_REF] Sparks | Potassium dynamics in soils[END_REF]. Recent evaluations of K-testing methods and their interpretation suggest however that common categories used to describe K in soils (soluble, exchangeable, fixed) are not universally agreed upon or understood and that renewed exploration of mineralogical controls over K cycling and interlayer K contributions to plant nutrition are needed [START_REF] Brouder | Potassium Cycling[END_REF].

Whilst exchangeable K, weakly held by the surface charges of clay minerals, and solution K are readily available for uptake by plants, release of mineral K which is present in the interlayer space of high-charge 2:1 clay minerals (micas) is thought to require depletion of exchangeable and solution K to low concentration [START_REF] Smith | Extractable Potassium in grundite illite: I. Method of extraction[END_REF]. Concentrations associated with such K release from high-charge 2:1 clay interlayers are known to vary with the initial mica crystal structure, ranging from 2.3 to 16.8 µg mL -1 for trioctahedral micas to <0.1 µg mL -1 for dioctahedral muscovite and illite [START_REF] Rausell-Colom | Studies in the artificial weathering of mica[END_REF][START_REF] Sparks | Potassium dynamics in soils[END_REF].

The former K concentrations are compatible with those reported in soil solutions (4-40 µg mL -1 , [START_REF] Hinsinger | How do plant roots acquire mineral nutrients? Chemical processes involved in the rhizosphere[END_REF] whereas the latter ones were only reported in the close vicinity (~1 mm) of plant roots where mineral dissolution may be significant [START_REF] Jungk | Ion Diffusion in the Soil-Root System[END_REF][START_REF] Jungk | Potassium depletion of the soil-root interface in relation to soil parameters and root properties[END_REF]. The possible release of mineral K + from mica interlayers and substitution by other cations is most often intimately related to clay mineralogical alteration, and more especially to the lowering of layer charge deficit [START_REF] Hinsinger | How do plant roots acquire mineral nutrients? Chemical processes involved in the rhizosphere[END_REF][START_REF] Marschner | Mineral nutrition of higher plants[END_REF][START_REF] Mortland | Alateraton of Biotite to Vermiculite by plant growth[END_REF][START_REF] Tributh | Effect of Potassium removal by crops on transformation of illitic clay minerals[END_REF][START_REF] Walker | Trioctahedral minerals in the soil-clays of North-East Scotland[END_REF]. Clay particle size may also influence the release of mineral K from clay minerals, which decreases with decreasing particle size [START_REF] Doll | Release of Potassium from Soil Fractions During Cropping[END_REF][START_REF] Reed | Kinetics of potassium release from biotite and muscovite in sodium tetraphenylboron solutions[END_REF]Scott and Reed, 1962b;[START_REF] Scott | Susceptibility of interlayer Potassium in micas to exchange with Sodium[END_REF].

The ability of plants to recover K from clay minerals has been observed in laboratory experiments from the vermiculitisation of K-bearing micas induced by plants grown on clay substrates [START_REF] Barré | Soil-plant potassium transfer: impact of plant activity on clay minerals as seen from X-ray diffraction[END_REF](Barré et al., , 2008a;;Hinsinger and Jaillard, 1993;[START_REF] Hinsinger | Rapid weathering of a trioctahedral mica by the roots of ryegrass[END_REF][START_REF] Niebes | Release of nonexchangeable potassium from different size fractions of two highly K-fertilized soils in the rhizosphere of rape (Brassica napus cv Drakkar)[END_REF]. The extent to which plants can extract K varies notably with plant species, proximity to plant roots, soil solution chemistry and clay mineralogy [START_REF] Barber | Root growth and phosphorus and potassium uptake by two corn genotypes in the field[END_REF][START_REF] Boyle | Biological weathering of silicate minerals[END_REF][START_REF] Kraffczyk | Soluble root exudates of maize: Influence of potassium supply and rhizosphere microorganisms[END_REF][START_REF] Marschner | Mineral nutrition of higher plants[END_REF][START_REF] Rengel | Crops and genotypes differ in efficiency of potassium uptake and use[END_REF][START_REF] Schenk | Potassium and phosphorus uptake by corn genotypes grown in the field as influenced by root characteristics[END_REF]Scott and Reed, 1962a). Long-term field experiments that began in the 19 th and 20 th centuries have led to the view that continuous cropping (30-135 years) increases smectite and/or vermiculite contents at the expense of mica/illite (Dissing Nielsen and [START_REF] Dissing Nielsen | The influence of K-depletion on mineralogical changes in pedons from two field experiments and in soils from four pot experiments[END_REF][START_REF] Møberg | Mineralogical changes in soils used for Potassium-depletion experiments for some years in pots and in the field[END_REF][START_REF] Tributh | Effect of Potassium removal by crops on transformation of illitic clay minerals[END_REF]. Surprisingly, [START_REF] Singh | Changes with time in the potassium content and phyllosilicates in the soil of the Broadbalk continuous wheat experiment at Rothamsted[END_REF] observed an increase in the average Kconcentration of soils from the Broadbalk continuous wheat plots established in 1843 at the Rothamsted Experimental Station (UK). The average total K-concentration increased from 1.01 % (in 1856) to 1.10 % and 1.20 % in 1987 in non K-fertilised and K-fertilised plots, respectively, but no significant modification of clay mineralogy was detected.

In the aforementioned studies, clay mineralogy was determined using qualitative Xray diffraction (XRD) methods that rely on the description of sharp and well-defined peaks associated with discrete clay minerals. Such qualitative description of XRD data does not recover comprehensive information about clay mineralogy and fails to measure contributions from mineral fractions potentially highly reactive: i) extremely small or cryptocrystalline and ii) complex interstratified minerals are essentially overlooked owing to their weakly modulated and faint diffraction signatures [START_REF] Lanson | Modelling of X-ray diffraction profiles: Investigation of defective lamellar structure crystal chemistry[END_REF]. Over the last two decades, the peak decomposition method has improved and rationalised interpretation of XRD data, by decomposing profiles into elementary contributions that allow identifying contributions from both discrete and interstratified clay minerals [START_REF] Lanson | Decomposition of experimental X-ray diffraction patterns (profile fitting): A convenient way to study clay minerals[END_REF]. Composition of the latter may be derived from the comparison with calculated patterns (e.g., [START_REF] Reynolds | NEWMOD: A computer program for the calculation of one-dimensionnal patterns of mixed-layered clays[END_REF]. As a result, shifts in the position of diffraction maxima and/or composition of mixed layers have been correlated to changes in K-concentration and to K-fertilisation (Barré et al., 2008a(Barré et al., , 2008b;;[START_REF] Pernes-Debuyser | Soil mineralogy evolution in the INRA 42 Plots experiment (Versailles, France)[END_REF][START_REF] Velde | Clay mineral changes in the Morrow experimental plots[END_REF]. In their study of the Morrow Plots long-term experiment, [START_REF] Velde | Clay mineral changes in the Morrow experimental plots[END_REF] reported with this approach the apparent illite-depletion of the clay fraction from the continuous corn plot, whereas no evolution was detectable for the three year rotation (corn-oats-hay) plot. As calculated from yearly yields and fertiliser inputs, both plots exhibit similar K balance over the recent period, however (Table 1 - [START_REF] Khan | The potassium paradox: Implications for soil fertility, crop production and human health[END_REF]. On the other hand, these authors reported mineral stability in both the unfertilised crop rotation plot and the fertilised continuous corn plots despite negative and positive K balances, respectively (Table 1 - [START_REF] Khan | The potassium paradox: Implications for soil fertility, crop production and human health[END_REF]. In any case, despite improving the description of clay mineral phases to include mixed layers and providing a proxy to their composition, the peak decomposition method does not account for the whole XRD profile and the obtained information remains essentially qualitative. In addition, the decomposition method does not allow taking into account highly-disordered and/or interstratified contributions with small particle size.

More recently, a full-profile modelling method has been developed and used to characterise clay mineralogy of hydrothermal and diagenetic sequences [START_REF] Aplin | Influence of mechanical compaction and clay mineral diagenesis on the microfabric and pore-scale properties of deep-water Gulf of Mexico Mudstones[END_REF][START_REF] Drits | Sequential structure transformation of illitesmectite-vermiculite during diagenesis of Upper Jurassic shales from the North Sea and Denmark[END_REF][START_REF] Inoue | Illite-smectite mixed-layer minerals in the hydrothermal alteration of volcanic rocks: I. One-dimensional XRD structure analysis and characterization of component layers[END_REF][START_REF] Lindgreen | Illite-smectite structural changes during metamorphism in black Cambrian Alum shales from the Baltic area[END_REF][START_REF] Mccarty | Heterogeneous mixed-layer clays from the Cretaceous greensand, Isle of Wight, southern England[END_REF]Sakharov et al., 1999aSakharov et al., , 1999b)). This approach has been adapted for soil samples and allows investigating soil mineralogy in a quantitative manner [START_REF] Hubert | Advances in characterization of soil clay mineralogy using X-ray diffraction: from decomposition to profile fitting[END_REF][START_REF] Hubert | Unraveling complex < 2 µm clay mineralogy from soils using X-ray diffraction profile modeling on particle-size sub-fractions: Implications for soil pedogenesis and reactivity[END_REF][START_REF] Viennet | Investigation of clay mineralogy in a temperate acidic soil of a forest using X-ray diffraction profile modeling: Beyond the HIS and HIV description[END_REF].

Quantitative mineralogy appears especially important in the perspective of expected increases in fertilisation and competition for arable land resources, but also for the growing interest in CO2 sequestration in cultivated soils [START_REF] Defries | Deforestation driven by urban population growth and agricultural trade in the twenty-first century[END_REF][START_REF] Lambin | Global land use change, economic globalization, and the looming land scarcity[END_REF][START_REF] Mueller | Closing yield gaps through nutrient and water management[END_REF]. It is evident that comprehensive and quantitative descriptions of soil (clay) mineralogy are key to the assessment of soil responses to land-use changes or fertiliser inputs, given the contradictory results obtained from studies based on qualitative analyses of XRD data (see for example Dissing Nielsen and [START_REF] Dissing Nielsen | The influence of K-depletion on mineralogical changes in pedons from two field experiments and in soils from four pot experiments[END_REF][START_REF] Singh | Changes with time in the potassium content and phyllosilicates in the soil of the Broadbalk continuous wheat experiment at Rothamsted[END_REF][START_REF] Tributh | Effect of Potassium removal by crops on transformation of illitic clay minerals[END_REF]. The aim of the present study was thus to obtain a comprehensive compositional description of the clay minerals from the Morrow Plots by using full-profile XRD modelling and to compare and contrast, quantitatively, the effects of 110 years of continuous plant-growth under different agronomic treatments on the nature and potassium content of clay minerals in the soil.

Materials and Methods

Location and characteristics

The Morrow Plots (MP) are located on the Urbana-Champaign campus of the University of Illinois (USA). The soil is an Aquic Argiudoll (loessic Flanagan silt loam -USDA soil taxonomy) developed on glacial till under temperate grassland vegetation. The MP receive an average of 968 mm precipitation (576 mm snowfall) per year, with a yearly average temperature of 11.1 °C (1904( -2014[START_REF] Mrcc | Cli-MATE: Midwestern Regional Climate Centre Application Tools Environment[END_REF]. The MP have been under continuous cultivation since their establishment in 1876 and are accompanied by an extensive record of crop species, rotations, fertiliser inputs and yields [START_REF] Aref | Long-term trends of corn yield and soil organic matter in different crop sequences and soil fertility treatments on the Morrow Plots[END_REF][START_REF] Odell | Changes in organic Carbon and Nitrogren of Morrow plots soils under different treatments, 1904-1973[END_REF]. Ten plots were established initially, and of these only three now remain. Each of the remaining plots is subject to a different crop rotation. Plot #3 is planted with corn (Zea mays) every year (C plots). Plot #4 (not included in this study) was initially planted in a two-yearly corn-oats (Avena sativa) rotation before being changed to a twoyearly corn-soybeans (Glycine max) rotation in 1967 to match regional cropping practices.

Plot #5 is planted in corn-oats-hay on a three-yearly rotation (R plots), where the hay crop was red clover (Trifolium pratense) from 1904 to 1954 and alfalfa (Medicago sativa) from 1955 onwards. Corn hybrids have been introduced in the rotations since 1937. Each of the plots is now divided into eight subplots, which each have a different fertilisation history. In the present study, samples were obtained from subplots which have remained unfertilised (CU and RU subplots) from the start of the experiment in 1876, and from subplots unfertilised between 1876 and 1954 which have received a nitrogen-phosphate-potassium (NPK) fertilisation treatment since 1955 (CF and RF subplots). Fertiliser is applied every year in CF subplot and once every three years (when corn is planted) in RF subplot. Fertilisation allowed for positive K balances in both CF and RF subplots, whereas [START_REF] Khan | The potassium paradox: Implications for soil fertility, crop production and human health[END_REF] computed negative balances for the equivalent unfertilised subplots (CU and RU subplots -Table 1). Even for unfertilised subplots, soil K tests exhibit rather high values for soils cropped continuously for over a century (Table 1 - [START_REF] Aref | Long-term trends of corn yield and soil organic matter in different crop sequences and soil fertility treatments on the Morrow Plots[END_REF]. The seeding density of F subplots was also increased in 1955. Only minor adjustments have taken place since 1955, most notably a change from mouldboard to chisel plough in 1989 and an increase in seeding density in 2013. In 2013 and 2014, sampling took place from the 0-25 cm surface horizon on fallow plots just after harvest (late-September to October). Unfertilised samples from 1904 Plots #3 and #5, and both fertilised and unfertilised subplots from 1957, 1980 and 1997 were collected from the soil archive at the U. of Illinois -Department of Agriculture (0-20 cm horizon). The latter samples were initially taken in the fall after harvest (Bob Dunker-Univ. of Illinois, personal communication). A complete description of the MP experimental fields can be found in [START_REF] Odell | Changes in organic Carbon and Nitrogren of Morrow plots soils under different treatments, 1904-1973[END_REF] and [START_REF] Aref | Long-term trends of corn yield and soil organic matter in different crop sequences and soil fertility treatments on the Morrow Plots[END_REF], where the samples in the current study are from NB (F) and NC (U) subplots of Plots #3 (C) and #5 (R).

Size-fractionation

Air dried samples were first sieved at 2000 µm to remove large organic fragments.

Approximately 5 g, accurately weighed, of <2000 µm soil material were dispersed via sonication in Milli-Q water (18.2 MΩ cm -1 ) before being wet-sieved to 50 µm (10 g of starting material were used for 2013 and 2014 samples). The 2000-50 µm fraction obtained was dried at 50 °C. The <50 µm fraction was centrifuged and flocculated with NaCl to reduce water content, then subjected to three overnight NaCl (2M) saturation cycles, followed by centrifugation and washing in dialysis tubes (Spectre/Por 6-8 kD membrane) until chloride free. Milli-Q water was systematically used for solution preparation, washing and dispersion steps of sample fractionation and preparation. Sequential size fractionation was then performed by repeated centrifugation using a Beckmann-Coulter Avanti J-20XP centrifuge (procedure adapted from [START_REF] Hubert | Unraveling complex < 2 µm clay mineralogy from soils using X-ray diffraction profile modeling on particle-size sub-fractions: Implications for soil pedogenesis and reactivity[END_REF]. A JS 4.3 swinging bucket rotor was used to separate 50-2 and <2 µm fractions, with the 50-2 µm fraction set aside to dry in the oven at 50 °C. 10 repetitions were performed to obtain a clear supernatant, with dispersion via sonication (30 s) between successive cycles. 2-0.2, 0.2-0.05 and <0.05 µm subfractions were obtained from repeated centrifugation of the initial <2 µm fraction using a JS 24.38 swinging bucket rotor, with 10-15 cycles being necessary to obtain a clear supernatant for each subfraction. 1 M CaCl2 or KCl was then used to saturate subfractions with either Ca 2+ or K + (three 24 h cycles) and samples were dialysed until chloride-free.

Chemical analyses

Bulk <2 µm fractions from 1904, 1957and 2014(CF and RU subplots -CU and RU for 1904) were analysed for major and trace elements by inductively coupled plasma optical emission after sample preparation according to [START_REF] Carignan | Routine Analyses of Trace Elements in Geological Samples using Flow Injection and Low Pressure On-Line Liquid Chromatography Coupled to ICP-MS: A Study of Geochemical Reference Materials BR, DR-N, UB-N, AN-G and GH[END_REF]. These samples were analysed also for C, H, and N to determine their organic composition. Cation exchange capacity (CEC) analysis was also performed on these bulk <2 µm fractions, and on 2-0.2, 0.2-0.05, and <0.05 µm subfractions of 2014 samples (CF and RU subplots) using the cobalt hexamine method [START_REF] Dohrmann | Three new, quick CEC methods for determining the amounts of exchangeable Calcium cations in calcareous clays[END_REF][START_REF] Orsini | Utilisation du chlorure de cobaltihexamine pour la détermination simultanée de la capacité d'échange et des bases échangeables des sols[END_REF].

X-ray diffraction analysis

X-ray diffraction (XRD) analysis of <2 µm subfractions was performed on Casaturated oriented slides prepared via the pipette method, under air-dried (Ca-AD) or ethylene glycol (Ca-EG) saturated conditions. The mass of sample present on each slide was determined and used for subsequent XRD profile modelling. EG solvation was ensured through vapour exposure at 50 °C overnight. Oriented slides of selected K-saturated samples were also analysed following heat treatments at 150 or 350 °C for four hours (K-150 and K-350 samples, respectively). Samples were allowed to cool to room temperature before XRD data collection. XRD patterns of Ca-AD and Ca-EG samples were collected at 40 % relative humidity using a Bruker D8 equipped with an Anton Paar CHC+ chamber coupled to an MHG Messtechnik humidity controller. Patterns of K-150 and K-350 samples were collected at 5 % relative humidity with the same instrument. A Sol-XE Si (Li) solid state detector from Baltic Scientific Instruments was used to record intensities using CuKα radiation over 2-50 °2Θ using 0.04 ° steps-and a counting time of 6 s per step for <2 µm subfractions. 50-2 μm fractions were prepared for randomly oriented powder XRD analysis by milling with ethanol in a McCrone® micronising mill for 8 minutes. XRD data was collected over 3-90 °2Θ with 0.026 ° steps and a 8 s counting time per step using the Bruker D8 diffractometer. In all cases, the sizes of the divergence slit, the two Soller slits and the antiscatter slit were 0.3 °, 2.5 °, 2.5 ° and 0.3 °, respectively.

XRD full-profile modelling

Qualitative analysis of the <2 µm subfractions was carried out as described in [START_REF] Bakker | Soil development under continuous agriculture at the Morrow Plots experimental fields from X-ray diffraction profile modelling[END_REF] where a detailed description of phase identification procedure is found.

Full-profile XRD pattern modelling was then performed with the Sybilla programme developed by Chevron™ [START_REF] Aplin | Influence of mechanical compaction and clay mineral diagenesis on the microfabric and pore-scale properties of deep-water Gulf of Mexico Mudstones[END_REF]. Parameterisation and fitting procedure described in [START_REF] Bakker | Soil development under continuous agriculture at the Morrow Plots experimental fields from X-ray diffraction profile modelling[END_REF] were used for all samples and subplots from 1904-2014.

Briefly, this modelling is based on the direct comparison between experimental and calculated 00l XRD reflections, and fitting takes place via trial-and-error determination of structure factors to arrive at a single structure model which fits both the Ca-AD and Ca-EG XRD profiles of <2 µm subfractions [START_REF] Drits | Sequential structure transformation of illitesmectite-vermiculite during diagenesis of Upper Jurassic shales from the North Sea and Denmark[END_REF][START_REF] Hubert | Advances in characterization of soil clay mineralogy using X-ray diffraction: from decomposition to profile fitting[END_REF][START_REF] Hubert | Unraveling complex < 2 µm clay mineralogy from soils using X-ray diffraction profile modeling on particle-size sub-fractions: Implications for soil pedogenesis and reactivity[END_REF]Sakharov et al., 1999aSakharov et al., , 1999b;;[START_REF] Viennet | Investigation of clay mineralogy in a temperate acidic soil of a forest using X-ray diffraction profile modeling: Beyond the HIS and HIV description[END_REF]. Such an approach relies on differences in layer-tolayer distances and structure factors for expandable layers arising from the different treatments (Ca-AD and Ca-EG) and leads to the determination of the proportion or weighted concentration of each independent phase in the sample.

Both discrete and interstratified contributions were considered. In the latter, factors such as layer types, their proportions and stacking sequences were refined. A single structure model which fits both the Ca-AD and Ca-EG XRD profiles of <2 µm subfractions is obtained at the end of the procedure. Discrete chlorite, kaolinite, illite and smectite were identified in <2 µm subfractions, in addition to up to seven mixed-layer contributions. The Profex interface for BGMN software [START_REF] Doebelin | Profex: a graphical user interface for the Rietveld refinement program BGMN[END_REF] was used for quantitative phase analysis through the Rietveld refinement of XRD patterns from the 50-2 µm fractions.

Results

Granulometry

Results of size-fractionation are presented in Fig. 1. Recovery was 90.7-97.8 % of the starting weight (94.9 % on average). Individual results show a high and non-systematic scatter between subplots for a given year, and for a single subplot in different years (Fig. S1).

As a result of this variability, plots were assessed on the basis of the cropping practices; yearly averages were also considered. When averaged over all subplots for a given year, the 2000-50 µm fraction represents 7 % (± 2 %), the 50-2 µm fraction 64 % (± 2 %), and the <2 µm fraction 29 % (± 2 %), for the years 1957 to 2014 (Fig. 1a). No systematic variation was observed between subplots for these fractions as a result of cropping practices (Figs. 1b-d, S1). These values are 4, 62, and 33 %, respectively, for 1904. The proportions of the 2-0.2, 0.2-0.05 and <0.05 µm subfractions are also significantly different for 1904 compared to other years, most likely indicating sample grinding, consistent with the absence of soil aggregates in the sampled jar. When looking at the relative proportions of <2 µm subfractions averaged over all subplots in a given year (Fig. 1d), there is a tendency for the <0.05 µm subfraction to increase from ~9 % of total sample mass in 1957 to ~13 % in 2013-2014 (Fig. 1d). At the same time the relative mass of the 0.2-0.05 µm subfraction decreases by ~2 % whereas that of the 2-0.2 µm subfractions is about stable (possible decrease by ~1 mass % -Figs. 1b-c). In addition, minor differences can be observed between the proportions of 2-0.2 µm subfractions of C and R subplots, C subplots containing generally slightly less mass than corresponding R subplots (Fig. 1b). The inverse is observed for <0.05 µm subfractions, where C subplots tend to have a greater proportion of <0.05 µm subfraction than corresponding R subplots, except for recent years. These differences are statistically insignificant over the period studied, however.

Chemical analyses

Elemental analyses (Table 2) are consistent with a soil clay fraction dominated by phyllosilicates, with high levels of SiO2 (43-46 %) and Al2O3 (17-19 %). K levels were similar in both subplots in 1904 (2.13 % K2O), and did not vary between 1904 and 1957 (2.15 %). These levels are similar to those reported for clay fractions from a variety of soils of similar origin [START_REF] Doll | Release of Potassium from Soil Fractions During Cropping[END_REF][START_REF] Singh | Changes with time in the potassium content and phyllosilicates in the soil of the Broadbalk continuous wheat experiment at Rothamsted[END_REF]. Whilst K decreased in clay fractions from both RU and CF subplots since 1957, the rate of decrease in RU subplots (-16 %) is faster than that in CF subplots (-7 %). In 1904, the C/N ratio of the RU subplot (11.3) was greater than that of CF (9.5). For both subplots, the C/N ratio has decreased since 1904 in line with the systematic loss of carbon. The carbon decrease rate has slowed in CF subplots since 1957 and the introduction of fertilisation, whilst C levels continue to decrease at a similar rate post-1957 in RU subplots. The CEC of the <2 µm fraction measured for 1904 is similar for CU and RU subplots and remains constant between 1904 and 1957 (Table 2). From 1957, the CEC of both C and R subplots increases significantly (by ~15 %) from 48 to 55 cmolc kg -1 in 2014 for both CF and RU subplots. In addition, CEC of the <0.05 µm subfractions are akin in these two subplots at 81-82 cmolc kg -1 (sample 2014) On the other hand, 0.2-0.05 and 2-0.2 µm subfractions demonstrate a difference of 7-8 cmolc kg -1 favouring RU subplots (Table 3).

XRD qualitative analysis and XRD profile simulation

Random powder XRD of 50-2 µm fractions and subsequent quantitative phase analysis found the mineralogy of the 50-2 µm fractions to be dominated by quartz (65-70 %) and potassium-(10-13 %) and plagioclase-feldspars (10-12 % -Fig. 2 and S2a-d, Table S1).

Mica, kaolinite and chlorite also contribute in minor amounts to this fraction (~ 5, 1 and 2 % respectively), whereas trace amounts (<1 %) of calcite, anatase and amphibole were also identified. There is no systematic variation in mineralogy of 50-2 µm fractions from subplotto-subplot or from year-to-year, with the exception of those of 1904 which had less mica and chlorite, presumably as the result of sample grinding (Table S1 and Figs. S2a-d). More especially, the relative abundances of K-bearing minerals (K-spars and micas) do not exhibit statistically significant evolutions with time and do not differ significantly as a function of agronomic practices (Table S1).

Full-profile modelling of the 2-0.2, 0.2-0.05 and <0.05 µm subfractions for CF, CU, RF and RU subplots from 1904, 1957[START_REF] Schenk | Potassium and phosphorus uptake by corn genotypes grown in the field as influenced by root characteristics[END_REF][START_REF] Lanson | Decomposition of experimental X-ray diffraction patterns (profile fitting): A convenient way to study clay minerals[END_REF], 2013and 2014 was performed using the structure model determined by [START_REF] Bakker | Soil development under continuous agriculture at the Morrow Plots experimental fields from X-ray diffraction profile modelling[END_REF] for the RU subplot. Similar contributions were used to model all <2 µm subfractions from 1904-2014, with only limited modifications of their composition (layer-type proportions) and coherent scattering domain (CSD) sizes (Tables S2 andS3). The <0.05 µm subfraction is dominated by randomly interstratified mixed layers in which illite, smectite and chlorite layers coexist (ISSCh). Two main contributions with contrasting illite contents (50 and 35 % in ISSCh 50 and ISSCh 35, respectively) coexist in this finest subfraction with two additional ISSCh contributions, one dominated by illite (ISSCh 90), the other by smectite (ISSCh 5 -Tables S2c andS3c). Minor amounts of discrete smectite with extremely small CSD sizes and of kaolinite-illite (KI) complement the structure model of this <0.05 µm subfraction. Except for discrete smectite, all contributions are present, with similar compositions and increased CSD sizes, in the 0.2-0.05 and 2-0.2 µm subfractions, together with discrete chlorite and discrete kaolinite. The peak at ~10 Å is accounted for by an additional illite-rich (ISSCh 80) mixed layer and discrete illite for the 0. S2 andS3. Experimental XRD patterns recorded for the same subfraction from different subplots or years show little variation in either peak position or intensity . Consistently, very little mineralogical difference was noted in the results of the full-profile modelling. Similarly, the composition of the different mixed layers varies only slightly between different subfractions. For example, for 2014 CF subplot ISSCh 50 has a 52/26/22 layer composition in the <0.05 µm subfraction which changes only slightly to 45/30/25 in the 2-0.2 µm subfraction (Ca-EG treatment -Table S3). This is fairly representative of the compositional evolution determined between fine and coarse clay subfractions for all subplots and phases. The content of smectite layers tends to decrease for 2-0.2 µm subfractions compared to finer subfractions whilst the proportion of illite layers decreases with increasing size subfraction for ISSCh 50, 35 and 5 but increases in ISSCh 90 and 80. These variations remain however marginal (~5 % of all layers at most). Similarly, composition of individual contributions appears stable as a function of time. For example, the composition of ISSCh 50 in the <0.05 µm subfractions varies from 55/23/22 (1904 CU), to 52/26/22 (1957 CF and 2014 CF), indicative of limited mineralogical evolution with time (Tables S2 andS3). Finally, K-spars are systematically present in the 2-0.2 µm clay subfraction, with no apparent evidence for an evolution of their relative contribution with time and/or agronomic practices (Figs. S3a-d).

Within each subfraction, proportions of the different layer types can be computed from the summation over all contributions of the proportion of this specific layer type within a contribution multiplied by the relative proportion of this contribution in the subfraction (Fig. 4). Proportions of illite, smectite and chlorite layers all increase in <0.05 µm subfractions from 1957 to 2014, whilst the amount of kaolinite is stable (Fig. 4a). There is a concurrent decrease in the proportions of illite, smectite and chlorite layers with time in 0.2-0.05 µm subfractions, kaolinite content being again essentially stable (Fig. 4b). The proportion of illite layers in 2-0.2 µm subfractions shows considerable scatter, making it difficult to define a clear trend. The proportion of illite layers is however stable or marginally decreasing with time, whilst smectite layers also show a marginal decrease and chlorite and kaolinite layers are stable. All evolutions observed remain too marginal given the time period studied (60 years) to be statistically significant, the uncertainty on mixed layer composition being about ±5% of a given layer type [START_REF] Hubert | Unraveling complex < 2 µm clay mineralogy from soils using X-ray diffraction profile modeling on particle-size sub-fractions: Implications for soil pedogenesis and reactivity[END_REF][START_REF] Lanson | Diagenetic Smectite-to-Illite Transition in Clay-Rich Sediments: a Reappraisal of X-Ray Diffraction Results Using the Multi-Specimen Method[END_REF]. These evolutions appear to be mostly linked to the evolution of gain size distribution within the <2 µm fraction (Fig. 1). When all clay subfractions are considered jointly (Fig. 4d) the bulk clay mineralogy appears unchanged in the MP over the last 60 years of cropping, despite some intrinsic scatter.

Discussion

Variability of results

The granulometry results show significant variability between subplots for a given year, this variation being non-systematic on a year-to-year-basis (Figs. 1a andS1). This variability is most likely related to heterogeneity within the MP, despite their small size, and is further indicative of the absence of systematic impact of crop rotation or fertilisation.

Sampling date scatter is possibly an additional source of variability that must be taken into consideration, as the exact timing of sampling post-harvest is not well-defined for historic samples (fallow period after harvest). Planting dates in the plots have varied by as much as 147 days during the MP experiment, suggesting that harvest and therefore sampling dates are also subject to considerable variation [START_REF] Aref | Long-term trends of corn yield and soil organic matter in different crop sequences and soil fertility treatments on the Morrow Plots[END_REF]. Given that numerous studies have related changes in soil properties over periods as short as two weeks [START_REF] Collignon | Seasonal dynamics of Al-and Fe-bearing secondary minerals in an acid forest soil: influence of Norway spruce roots (Picea abies (L.) Karst.)[END_REF][START_REF] Haines | Seasonal variation in properties of five forest soils in Southwest Georgia[END_REF][START_REF] Marschner | Mineral nutrition of higher plants[END_REF]Mengel et al., 2001;[START_REF] Peoples | Soil: Nutrient Cycling[END_REF][START_REF] Turpault | Clay minerals: Precise markers of the spatial and temporal variability of the biogeochemical soil environment[END_REF][START_REF] Weaver | Seasonal variation in soil nutrients under six Rocky Mountain vegetation types[END_REF], sampling date scatter may play a role in the observed granulometric and mineralogical variability. Furthermore, analysis of yield and climate data of the MP from 1888 to 1996 [START_REF] Aref | Long-term trends of corn yield and soil organic matter in different crop sequences and soil fertility treatments on the Morrow Plots[END_REF] revealed that yield was positively correlated not only with the temperature and precipitation of the preceding November and December but also with those of the preceding year. Climatic fluctuations also affect the distribution of plant roots in soils, and therefore the location of nutrient uptake, and this variation may be significant on a year-to-year basis [START_REF] Marschner | Mineral nutrition of higher plants[END_REF]Mengel et al., 2001). The impacts of such climatic fluctuations on soil moisture and K availability may account, at least partially, for the difference observed between 2013 and 2014 samples despite their small separation in time. Fertilisation treatments represent another source of variability. Whilst NPK fertilisation positively influenced average yields in F subplots, yield stability dropped when rates of fertiliser application were increased [START_REF] Aref | Long-term trends of corn yield and soil organic matter in different crop sequences and soil fertility treatments on the Morrow Plots[END_REF], reflecting variations in the short term availability of plant nutrients.

Evolution of (clay) mineralogy

Following ~150 years of continuous cultivation, mineralogical differences in the MP are limited to a statistically insignificant increase of illite layers in 2-0.2 µm subfractions from U subplots compared to F subplots and of smectite layers in <0.05 µm subfractions from C subplots compared to R ones (Fig. 4). In addition, K-spars are systematically present in the 2-0.2 µm subfraction, with no evidence of dissolution with time even in U subplots which exhibit strongly negative K-balances (Table 1; Figs. S3a-d), whereas the mineralogy of the coarser 50-2 µm is essentially independent of time and agronomic practices (Table S1 and Figs. S2a-d). Despite the presence of K-spars in clay subfractions, which is favourable to their dissolution [START_REF] Manning | Testing the ability of plants to access potassium from framework silicate minerals[END_REF][START_REF] Sadusky | Kinetics and Mechanisms of Potassium Release from Sandy Middle Atlantic Coastal Plain Soils1[END_REF], their contribution to K supply is thus likely limited. This long-term mineralogical stability is consistent with previous reports demonstrating that mineralogical differences reported between rhizosphere and bulk soil after 40 days of corn cultivation did not persist longer than a growth cycle (130 days - [START_REF] Adamo | Short term clay mineral release and re-capture of potassium in a Zea mays field experiment[END_REF]. This long-term resilience most likely indicates that soil K reserves exceed the plough layer and that a significant proportion of plant nutrition is obtained from subsoil in a non K-limited context [START_REF] Hinsinger | Root traits and rhizosphere characteristics determining Potassium acquisition from soils[END_REF][START_REF] Kautz | Nutrient acquisition from arable subsoils in temperate climates: A review[END_REF][START_REF] Khan | The potassium paradox: Implications for soil fertility, crop production and human health[END_REF][START_REF] Kuhlmann | Importance of the subsoil for the K nutrition of crops[END_REF][START_REF] Marschner | Mineral nutrition of higher plants[END_REF]. The absence of K-limitation is substantiated by the observed yield increase in U subplots following the introduction of hybrid crops in 1937 [START_REF] Aref | Long-term trends of corn yield and soil organic matter in different crop sequences and soil fertility treatments on the Morrow Plots[END_REF], increase allowed by the high K-levels measured in non-amended soils cropped continuously for decades (217 and 242 kg K ha -1 measured in 1955 in R and C subplots, respectively -Table 1). Furthermore, these soil K-test values systematically increase with time whatever the K-balance in the considered plot (Table 1). In addition to the likely subsoil contribution, the similarity between F and U subplots may be favoured by the higher seeding density in the former -20,000 and 30,000 plants ha -1 in CU and RU subplots, respectively, versus 40,000 and 60,000 plants ha -1 (for 1957-1977 and 1978-2013, respectively) in F subplots [START_REF] Odell | Changes in organic Carbon and Nitrogren of Morrow plots soils under different treatments, 1904-1973[END_REF].

The lack of mineralogical evolution with time reported by [START_REF] Bakker | Soil development under continuous agriculture at the Morrow Plots experimental fields from X-ray diffraction profile modelling[END_REF] for RU subplots is upheld in the present study that considers different subplots and treatments, thus indicating a minimal influence of agronomic practices on soil mineralogy. These results are in agreement with [START_REF] Velde | Clay mineral changes in the Morrow experimental plots[END_REF] for R subplots, but the latter authors reported a decrease in illite in the clay fraction of C subplots, based on XRD profile decomposition.

The subtle variations in the position and shape of the 001 reflections in XRD profiles that were tentatively quantified by [START_REF] Velde | Clay mineral changes in the Morrow experimental plots[END_REF] were also observed in the present study. These were not systematic however and did not resolve into significant mineralogical differences following full-profile modelling, but rather contributed to the intrinsic sample variability.

Overall, the main evolution of clay mineralogy with time is the increased proportion of the <0.05 µm subfraction from 1957 to 2014 for all subplots regardless of agronomic practices. This increase is accompanied by an increase of the bulk clay CEC owing to the major contribution of the finest subfraction to the CEC of the bulk <2 µm fraction (Table 3).

Contribution of soil organic matter to the CEC increase is likely minor despite decreasing C/N ratios with time. The minor character of this contribution is supported by similar CEC values measured for both CF and RU subplots despite the contrasting evolutions of their C/N ratios.

The observed granulometric evolution of the bulk <2 µm fraction is most likely due to mineral dissolution to supplement K for plant nutrition, as the K concentration decreases with time in this fraction for all subplots (Table 2). This decrease of K concentration is consistent with the negative K balances reported for U subplots over the period 1955[START_REF] Inoue | Illite-smectite mixed-layer minerals in the hydrothermal alteration of volcanic rocks: I. One-dimensional XRD structure analysis and characterization of component layers[END_REF][START_REF] Khan | The potassium paradox: Implications for soil fertility, crop production and human health[END_REF]. These authors calculated positive K balances for F subplots, however, whereas a decrease in total K is also observed, a seemingly contradictory result. The mineral contributions to all subfractions being essentially similar in composition (Tables S2 andS3), the observed evolution in grain size likely takes place through dissolution of the phases present in the coarser 2-0.2 and 0.2-0.05 µm subfractions rather than through their alteration or transformation. The breakdown of micro-aggregates initially present in these coarser subfractions could also account, at least partially, for the granulometric and CEC evolutions observed with time. This effect is likely minor, however, owing to the systematic presence of minor amounts of carbonates in all soil samples despite the slightly acidic soil pH (Table S1). Furthermore, the steady decrease in C concentration from 1904-2014 is essentially disconnected from the CEC increase, that occurs mainly from 1957-2014 (Table 2). Protection of <0.05 µm particles by soil organic matter and sesquioxides, which are preferentially associated with this subfraction (e.g., [START_REF] Barré | Impact of phyllosilicate mineralogy on organic carbon stabilization in soils: incomplete knowledge and exciting prospects[END_REF][START_REF] Lützow | Stabilization of organic matter in temperate soils: mechanisms and their relevance under different soil conditions -a review[END_REF][START_REF] Singh | Stabilization of soil organic Carbon as influenced by clay mineralogy[END_REF], likely enhances the stability of this finest subfraction.

Mechanism of mineral K uptake by plants: Exchange Vs. dissolution

As discussed above, both the absence of significant modification of clay mineralogy and the increased abundance of the finest clay subfraction (<0.05 µm) at the expense of the coarser ones (0.2-0.05 and 2-0.2 µm subfractions) suggest that in the MP experiment K is most likely released from illite interlayer through dissolution rather than by exchange process. This inference raises the question of the prevalence of one process over the other and, more specifically, of the conditions leading to this predominance. Diffusion of K in illite interlayers at equilibrium in solutions similar to natural ones [START_REF] De Haan | Diffusion of Potassium-40 into an illite during prolonged shaking[END_REF][START_REF] Nye | The measurement and mechanism of ion diffusion in soils. VIII -A theory for the propagation of changes of pH in soils[END_REF] implies a negligible contribution of interlayer K exchange. However, early studies have unambiguously shown that K can be removed from K-rich mica interlayers (muscovite, illite, phlogopite, biotite) through the sole exchange process [START_REF] Reed | Kinetics of potassium release from biotite and muscovite in sodium tetraphenylboron solutions[END_REF][START_REF] Scott | Mechanisms of potassium release by soil minerals[END_REF][START_REF] Scott | Mechanisms of potassium release by soil minerals[END_REF]Reed, 1962a, 1962b;[START_REF] Scott | Susceptibility of interlayer Potassium in micas to exchange with Sodium[END_REF][START_REF] Smith | Extractable Potassium in grundite illite: I. Method of extraction[END_REF]. In this case, the driving force is the chemical gradient between mica interlayers and the solution, and the limiting factor is thus the K concentration in solution; threshold concentrations in solution allowing K exchange range from 2.3 to 16.8 µg mL -1 for trioctahedral micas to <0.1 µg mL -1 for dioctahedral muscovite and illite [START_REF] Rausell-Colom | Studies in the artificial weathering of mica[END_REF][START_REF] Sparks | Potassium dynamics in soils[END_REF]. The higher content of structural Fe(II) in natural trioctahedral micas, compared to dioctahedral ones, may be responsible for these contrasting behaviours, as Fe(II) oxidation decreases layer charge deficit, thus favouring the exchange of interlayer K + cations. Although low, these values allow for K exchange from trioctahedral micas in the laboratory (Hinsinger and Jaillard, 1993;[START_REF] Hinsinger | Rapid weathering of a trioctahedral mica by the roots of ryegrass[END_REF][START_REF] Niebes | Release of nonexchangeable potassium from different size fractions of two highly K-fertilized soils in the rhizosphere of rape (Brassica napus cv Drakkar)[END_REF][START_REF] Norouzi | Ability of alfalfa (Medicago sativa L.) to take up potassium from different micaceous minerals and consequent vermiculitization[END_REF], and in the field (Kapoor, 1972b;[START_REF] Wilson | The weathering of biotite in some Aberdeenshire soils[END_REF]. The lower threshold concentrations determined for dioctahedral micas most likely hamper K exchange, however, consistent with the conclusions of Kapoor (1972aKapoor ( , 1972b)), although root activity may induce such low concentrations in the rhizosphere [START_REF] Hinsinger | Acquisition of phosphorus and other poorly mobile nutrients by roots. Where do plant nutrition models fail?[END_REF][START_REF] Hinsinger | Root traits and rhizosphere characteristics determining Potassium acquisition from soils[END_REF][START_REF] Jungk | Ion Diffusion in the Soil-Root System[END_REF][START_REF] Walker | Absorption of potassium and rubidium from the soil by corn roots[END_REF]. The concentrations of other nutrients (e.g., P and NH4 + ) in solution may also have a key influence on the exchange of interlayer K, which can be totally inhibited in the presence of ~5 µmol L -1 of NH4 + , a concentration common in natural soils [START_REF] Springob | Blocking the release of potassium from clay interlayers by small concentrations of NH4+ and Cs+[END_REF].

Depending on the aforementioned experimental conditions, the kinetic competition between dissolution and exchange processes accounts of the heterogeneity of the results obtained from laboratory and field studies. Hinsinger and coworkers [START_REF] Hinsinger | Rapid weathering of a trioctahedral mica by the roots of ryegrass[END_REF](Hinsinger et al., , 1993;;Hinsinger and Jaillard, 1993) showed unambiguously the prevalence of exchange leading to the fast vermiculitisation of coarse trioctahedral mica (2-100 µm phlogopite) by ryegrass and rape, possibly as the result of structural Fe(II) partial oxidation. Similar results were obtained by [START_REF] Naderizadeh | Organic matter induced mineralogical changes in clay-sized phlogopite and muscovite in alfalfa rhizosphere[END_REF] and [START_REF] Norouzi | Ability of alfalfa (Medicago sativa L.) to take up potassium from different micaceous minerals and consequent vermiculitization[END_REF] also for trioctahedral micas (phlogopite and biotite) under alfalfa growth, whereas [START_REF] Vetterlein | Illite transformation and potassium release upon changes in composition of the rhizophere soil solution[END_REF] reported a minor evolution of a reference illite (IMt-2, Clay Minerals Society Repository -<20 µm particles) to illlite-expandable interstratification following 100 days of alfalfa growth. Conversely, no mineralogical evolution was reported for a natural soil [START_REF] Vetterlein | Illite transformation and potassium release upon changes in composition of the rhizophere soil solution[END_REF] or muscovite [START_REF] Naderizadeh | Organic matter induced mineralogical changes in clay-sized phlogopite and muscovite in alfalfa rhizosphere[END_REF][START_REF] Norouzi | Ability of alfalfa (Medicago sativa L.) to take up potassium from different micaceous minerals and consequent vermiculitization[END_REF], under similar experimental conditions. Muscovite dissolution even appears to be limited as alfalfa displays K deficiency symptoms [START_REF] Norouzi | Ability of alfalfa (Medicago sativa L.) to take up potassium from different micaceous minerals and consequent vermiculitization[END_REF]. Prevalence of dissolution was observed by [START_REF] Feigenbaum | Dissolution of Illite-A possible mechanism of Potassium release[END_REF] using Ca and H-resins to release K from illite interlayers. Kinetic competition between dissolution and exchange processes was investigated directly by [START_REF] Viennet | Effect of particle size on the experimental dissolution and auto-aluminization processes of Kvermiculite[END_REF] in the specific context of K-vermiculite aluminisation under acidic conditions (trioctahedral high-charge clay mineral, pH = 3). These authors concluded that despite the faster release of both exchangeable (at the particle surface) and nonexchangeable interlayer (i.e. mineral) K for the finest size fraction investigated (0.1-0.2 µm), the exchange process was rapidly stopped. Release of interlayer K subsequently proceeds through stoichiometric dissolution of this finest size fraction, whereas exchange remains active for coarser fractions (1-2 µm and even more so for 10-20 µm). The decrease of K exchange when decreasing size fraction to clay size is consistent with experiments involving K release by chemical exchange [START_REF] Doll | Release of Potassium from Soil Fractions During Cropping[END_REF][START_REF] Reed | Kinetics of potassium release from biotite and muscovite in sodium tetraphenylboron solutions[END_REF]Scott and Reed, 1962b;[START_REF] Scott | Susceptibility of interlayer Potassium in micas to exchange with Sodium[END_REF] or by cropping [START_REF] Doll | Release of Potassium from Soil Fractions During Cropping[END_REF], whereas the increased dissolution rate of finer subfractions is supported by experimental studies [START_REF] Acker | The influence of pH on biotite dissolution and alteration kinetics at low temperature[END_REF][START_REF] Bray | The effect of pH, grain size, and organic ligands on biotite weathering rates[END_REF]. Results of the present study support a marginal release of interlayer K by exchange mechanisms and the prevalence of dissolution over exchange, consistent with this general discussion, as MP surface soils are typified by fine-grained dioctahedral micas and clay minerals (inset Fig. 2), both features being unfavourable to K exchange from mica interlayers.

Conclusions

Mineralogy of the silt fraction and of clay subfractions was determined quantitatively by modelling X-ray diffraction (XRD) patterns for a temporal series of top soil samples collected from the Morrow Plots long-term experimental fields. Samples from different subplots encompass a variety of cropping practices resulting in positive or negative K balances (fertilised and unfertilised subplots, respectively). Mineralogy of the different size fractions is essentially stable over time, independent of the overall K-balance and of agronomic practices, showing no significant decrease of the mica or K-spar relative abundances as the result of K uptake by plants. Within the clay fraction, the relative proportion of the finest (<0.05 µm) subfraction increases with time, together with the cation exchange capacity of the bulk <2 µm fraction. This systematic granulometric evolution results from the dissolution of illite-rich coarse clay fractions. Prevalence of dissolution over K-exchange to supply plants with needed K most likely results from the dioctahedral character of fine-grained micas present in the Morrow Plots. Subsoil provide additional K for plant nutrition, whereas contribution of K-feldspars appears minimal in a non K-limited context. Finally, the improved description of (clay) mineralogy allowed by the full-profile modelling approach did not allow to evidence significant variation of (clay) mineralogy possibly resulting from cropping practices and related K-uptake.

Figure Captions

  2-0.05 µm and the 2-0.2 µm subfractions. Examples of the fit for 2014 CF and RU samples are shown in Fig. 3 whereas complete results of full-profile modelling of clay subfractions are reported in Tables

Fig. 1 .

 1 Fig. 1. Grain size distributions obtained for 1904-2014 soil samples with the sequential

Fig. 2 .

 2 Fig. 2. Selection of quantitative phase analysis performed with the Rietveld method on XRD patterns of sample 2014 (CF and RU subplots). Experimental data is represented by grey crosses, and solid black lines indicate the calculated intensity. Difference plots are solid grey lines underneath. Inset b) shows the position of the 06,33 reflection for the bulk <2 µm fraction of sample 2014 (RU subplot). Colour figure is available online.

Fig. 3 .

 3 Fig. 3. Selection of XRD full-profile modelling of clay subfractions. a) 2014 RU (left) and b) 2014 CF (right) subplots. Ca-EG and Ca-AD are shown at the top and the bottom of the Figure, respectively. XRD patterns of 2-0.2, 0.2-0.05 and <0.05 µm subfractions are shown from top to bottom for each treatment. Symbols as in Fig. 2. Quartz, feldspar, and hk contributions from phyllosilicates are indicated by solid black, open black and solid grey triangles respectively. These contributions were not taken into account during the modelling procedure. The vertical grey bars represents an increase in intensity of the high-angle region compared to the 2-22° 2θ angular range (×2 scale factor). Colour figure is available online.

Fig. 4 .

 4 Fig. 4. Relative abundances of illite, smectite, chlorite and kaolinite layers in the different clay subfractions (a: <0.05 µm, b: 0.2-0.05 µm, and c: 2-0.2 µm, c: 0.2-0.05 µm) and in the bulk <2 µm fraction (d). Yearly averages are represented by lines, symbols representing different agronomic practices for illite. Smectite is represented in grey to enhance the difference between illite and smectite. Note the y-axis scale in d) is twice that of the other

  al. graphical abstract

Table 1

 1

Table 3

 3 Measured and calculated CEC values for sample 2014 clay-size fractions and subfractions. The calculated CEC value for the bulk <2 µm fraction is a weighted average of the CEC values measured for the different subfractions.

			2014 RU			2014 CF
	Subfraction	CEC (meas.) cmolc kg -1	Rel. prop. (wt %)	CEC (calc.) cmolc kg -1	CEC (meas.) cmolc kg -1	Rel. prop. (wt %)	CEC (calc.) cmolc kg -1
	<2 µm	55.2	-	57.8	55.3	-	52.9
	2-0.2 µm	22.1	30.7	-	14.0	31.8	-
	0.2-0.05 µm	58.6	25.6	-	51.3	23.0	-
	<0.05 µm	82.2	43.8	-	81.2	45.1	-

Acknowledgments

Martine Lanson (ISTerre, Grenoble, France) has coordinated all chemical analyses and has performed CEC measurements. Chemical analyses were performed at the Service d'Analyse des Roches et des Minéraux (SARM-CNRS, Vandoeuvre-lès-Nancy, France) for major elements and at the Institut de Chimie des Substances Naturelles (ICSN -Gif sur Yvette, France) for C, H and N. Paolo Benavides (ISTerre, Grenoble, France) is thanked for size fractionation and XRD analysis of 1980 and 1997 samples. Tauhid B. Khan (ISTerre, Grenoble, France) is thanked for preliminary investigations of MP samples. Financial support from INSU/EC2CO program to FH and BL is acknowledged (Project Claie, 2015-2016). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56).

CU

CF RU RF

Exchangeable K content in surface soil (0-15 cm) from 1955 (kg K ha -1 ) a 242 217

Soil K content