Tchodjowiè P I Kpemoua

Sarah Leclerc

Pierre Barré

Sabine Houot

Valérie Pouteau

Cédric Plessis

Claire Chenu
email: claire.chenu@inrae.fr

Are carbon-storing soils more sensitive to climate change? A laboratory evaluation for agricultural temperate soils

Keywords: Carbon mineralization, agroecological practices, organic waste products, climate change, sensitivity

Soil moisture regime and temperature significantly affect the carbon mineralization • Dry wet cycles did not stimulate the carbon mineralization relative to wet controls

• Specific carbon mineralization did not differ among the three cropping systems at La Cage • Lower specific carbon mineralization in soils receiving organic waste products at QualiAgro

• Carbon-storing soils have a similar sensitivity to climate events as baseline soils

Introduction

The 4p1000 initiative encourages the establishment of agricultural practices that increase and/or preserve soil carbon stocks (www.4p1000.org). At the field scale, changes in soil organic carbon (SOC) stocks result from an imbalance between carbon inputs (crop residues, litterfall, root exudates, exogeneous organic matter application, etc.) and carbon outputs due to SOC mineralization, leaching or erosion [START_REF] Lal | Digging deeper: A holistic perspective of factors affecting soil organic carbon sequestration in agroecosystems[END_REF]. It is generally agreed that the most efficient way to increase SOC stocks is to increase carbon inputs (e.g., [START_REF] Virto | Carbon input differences as the main factor explaining the variability in soil organic C storage in no-tilled compared to inversion tilled agrosystems[END_REF][START_REF] Autret | Alternative arable cropping systems: A key to increase soil organic carbon storage? Results from a 16 year field experiment[END_REF][START_REF] Fujisaki | Data synthesis of carbon distribution in particle size fractions of tropical soils: Implications for soil carbon storage potential in croplands[END_REF]. This can be achieved by increasing field biomass production and residue return (e.g., cover crops, [START_REF] Poeplau | Carbon sequestration in agricultural soils via cultivation of cover crops -A meta-analysis[END_REF][START_REF] Autret | Alternative arable cropping systems: A key to increase soil organic carbon storage? Results from a 16 year field experiment[END_REF], or by mobilizing and spreading external carbon resources such as organic waste products (OWPs) [START_REF] Peltre | RothC simulation of carbon accumulation in soil after repeated application of widely different organic amendments[END_REF][START_REF] Paetsch | Urban waste composts enhance OC and N stocks after long-term amendment but do not alter organic matter composition[END_REF]. The implementation of agroecological practices that increase carbon inputs and/or reduce carbon outputs (e.g., conservation agriculture, organic agriculture, OWPs application) generally result to additional carbon storage [START_REF] Peltre | RothC simulation of carbon accumulation in soil after repeated application of widely different organic amendments[END_REF][START_REF] Autret | Alternative arable cropping systems: A key to increase soil organic carbon storage? Results from a 16 year field experiment[END_REF][START_REF] Paetsch | Urban waste composts enhance OC and N stocks after long-term amendment but do not alter organic matter composition[END_REF][START_REF] Pellerin | Stocker du carbone dans les sols français[END_REF]; such soils are referred to "carbon storing soils" in this study. However, the question is whether the adoption of these agroecological practices would still be beneficial for carbon sequestration in a warmer, drier, wetter climate or with intense drywet cycles?

It is widely recognized that climate change affects carbon mineralization in soils [START_REF] Joly | Resolving the intricate role of climate in litter decomposition[END_REF]. Climate change is leading to an overall temperature increase, as well as increases in the frequency and intensity of extreme climatic events such as droughts, heavy precipitation, fires, or freeze-thaw cycles [START_REF] Shukla | Climate Change and Land: an IPCC special report on climate change, desertification, land degradation, sustainable land management, food security, and greenhouse gas fluxes in terrestrial ecosystems[END_REF]. There is also an increase in the evaporation of terrestrial water and an increase in the moisture level in the atmosphere, leading to frequent cycles of drought and heavy rainfall events [START_REF] Dai | Increasing drought under global warming in observations and models[END_REF][START_REF] Donat | More extreme precipitation in the world's dry and wet regions[END_REF]. Overall, these climatic events could have positive or negative impacts on CO2 emissions into the atmosphere by affecting ecological processes which control the dynamics of organic carbon in agricultural soils. The existence of this important feedback effect has prompted several studies of how climate change affects the kinetics of SOC mineralization and the emission of CO2 from soils [START_REF] Kirschbaum | Will changes in soil organic carbon act as a positive or negative feedback on global warming?[END_REF][START_REF] Smith | Greenhouse gas mitigation in agriculture[END_REF]. For example, increased soil temperature is known to accelerate SOC decomposition because temperature-dependent reactions performed by microorganisms [START_REF] Davidson | Temperature sensitivity of soil carbon decomposition and feedbacks to climate change[END_REF] result in more rapid CO2 emissions from soil to the atmosphere (Trumbore andCzimczik, 2008, Karhu et al., 2014). Soil moisture can have a great impact on SOC decomposition by affecting the oxygen diffusion into the soil and the substrate availability for soil microorganisms [START_REF] Linn | Effect of Water-Filled Pore Space on Carbon Dioxide and Nitrous Oxide Production in Tilled and Nontilled Soils[END_REF][START_REF] Suseela | Effects of soil moisture on the temperature sensitivity of heterotrophic respiration vary seasonally in an old-field climate change experiment[END_REF], Moyano et al., 2013, Wang et al., 2014[START_REF] Zhou | Effects of Soil Moisture on the Temperature Sensitivity of Soil Heterotrophic Respiration: A Laboratory Incubation Study[END_REF][START_REF] Sierra | Sensitivity of decomposition rates of soil organic matter with respect to simultaneous changes in temperature and moisture[END_REF]. With climate change, it is expected that drywet cycles will increasingly occur across the globe. Droughts would tend to decrease biological activity, but the subsequent precipitation would lead to "pulses" of CO2 emissions to the atmosphere [START_REF] Birch | The effect of soil drying on humus decomposition and nitrogen availability[END_REF]. Understanding the effects of temperature and soil moisture regime on SOC dynamics, especially in carbon-storing soils, and assessing the sensitivity of carbon mineralization to these climatic events is necessary to anticipate future carbon losses in terrestrial ecosystems under climate change. Climate change will obviously also play on C inputs to soil which will influence SOC stock evolution, but this is outside the scope of our present study. [START_REF] Conant | Temperature and soil organic matter decomposition rates -synthesis of current knowledge and a way forward[END_REF] define temperature sensitivity as the rate of a process (decomposition, desorption) at a given temperature compared to a control temperature. Temperature sensitivity, determined by Q10, represents a proportional increase (or decrease) in SOC mineralization for a temperature increase of a 10 °C [START_REF] Kirschbaum | The temperature dependence of soil organic matter decomposition, and the effect of global warming on soil organic C storage[END_REF]. The Q10 is an important parameter for predicting the fates of SOC under warmer climate [START_REF] Kan | Carbon mineralization and its temperature sensitivity under no-till and straw returning in a wheat-maize cropping system[END_REF]. Findings of other studies showed that Q10 values varied largely across the range of temperatures [START_REF] Ghosh | Temperature sensitivity of soil organic carbon decomposition as affected by long-term fertilization under a soybean based cropping system in a sub-tropical Alfisol[END_REF][START_REF] Ghosh | Long-term fertilization effects on 13 C natural abundance, soil aggregation, and deep soil organic carbon sequestration in an Alfisol[END_REF]. Contrary to temperature, there is no standardized parameter to assess the sensitivity of carbon mineralization to soil moisture regime. In this study, we defined moisture sensitivity as the proportional response of soil microbial respiration (SOC mineralization) of soils exposed to different moisture levels relative to a control moisture. In addition to climatic factors, the quality of soil organic matter should also affect the rate of decomposition [START_REF] Guntiñas | Sensitivity of soil respiration to moisture and temperature[END_REF]. In soils formed under similar climates and from similar parent material, the quality of the organic matter depends on the vegetation cover, the quality of external organic input and the type of use and management to which the soil is subjected. So, the implementation of agroecological practices including conservation agriculture, organic agriculture and organic waste products application (composts, manure) affects the quantity and quality of organic matter in soils which, in turn, can play on the sensitivity of SOC mineralization to climate change.

No-tillage and crop rotation types under conservation agriculture affect the distribution of carbon in soil fractions (Zhang et al., 2020b). Some research suggests that additional C from increased residue inputs accumulates mostly in particulate organic matter (POM) fractions that are on average easily mineralized and only small gains occur in the resistant C pool [START_REF] Bhattacharyya | Long term fertilization effects on soil organic carbon pools in a sandy loam soil of the Indian sub-Himalayas[END_REF][START_REF] Mitchell | Amount and incorporation of plant residue inputs modify residue stabilisation dynamics in soil organic matter fractions[END_REF]. Long-term organic waste supply significantly increased the POM fractions [START_REF] Peltre | RothC simulation of carbon accumulation in soil after repeated application of widely different organic amendments[END_REF][START_REF] Paetsch | Urban waste composts enhance OC and N stocks after long-term amendment but do not alter organic matter composition[END_REF]. These results show that the implementation of agroecological practices increasing C inputs may affect SOC temporal stability. However, there is a lack of consensus in the literature on how SOC accumulated under alternative practices respond to seasonal climate variability and projected climate change (Davidson andJanssens, 2006, Carey et al., 2016). The carbon loss through heterotrophic respiration is intrinsically related to the availability and quality of this organic carbon [START_REF] Lindén | Contrasting effects of increased carbon input on boreal SOM decomposition with and without presence of living root system of Pinus sylvestris L[END_REF], Hopkins et al., 2014). The POM fraction is believed to be a quite labile carbon pool in soils, with a mean residence time up to 20 years [START_REF] Balesdent | Treize ans d'investigation de la dynamique du carbone organique des sols au moyen des abondances naturelles en 13C : bilan et perspectives[END_REF]. Therefore, if additional SOC carbon mainly comprises POM, as suggested in the literature, then the sensitivity of SOC mineralization to temperature and moisture regime under agroecological practices (carbon-storing soils) could differ from that of soils where agroecological practices have not been implemented (baseline soils). Indeed, previous studies have shown that labile carbon is less sensitive to temperature than stable carbon [START_REF] Liski | CO₂ Emissions from soil in response to climatic warming are overestimated: The decomposition of old soil organic matter is tolerant of temperature[END_REF][START_REF] Conen | Temperature sensitivity of young and old soil carbon -Same soil, slight differences in 13C natural abundance method, inconsistent results[END_REF][START_REF] Lefèvre | Higher temperature sensitivity for stable than for labile soil organic carbon -Evidence from incubations of long-term bare fallow soils[END_REF][START_REF] Xu | Distinct temperature sensitivity of soil carbon decomposition in forest organic layer and mineral soil[END_REF]. Based on these findings, we hypothesized that SOC mineralization in carbon-storing soils would be relatively less sensitive to temperature. Conversely, soil moisture affects the oxygenation of these POMs and their accessibility to decomposers by increasing the probability of contact between organic substrates and microorganisms [START_REF] Monard | Impact of soil matric potential on the fine-scale spatial distribution and activity of specific microbial degrader communities[END_REF]. Thus, the sensitivity of SOC mineralization to soil moisture regime is expected to be relatively higher in carbon-storing soils than in their baseline soils. Testing these hypotheses is crucial to provide knowledge that will allow to identify the best management practices that reduce CO2 emissions and enhance the soil carbon sink effect, even under climate change.

We aimed to investigate in this study two related research questions. The first objective was to assess the sensitivity of SOC mineralization to the soil moisture regime in carbon-storing soils compared to baseline soils. To achieve this objective, we incubated in the laboratory soils collected from longterm field experiments at different matric potentials (pF 1.5; pF 2.5 and pF 4.2) and subsequent dry (pF 4.2)-wet (pF 1.5) cycles at 20°C and monitored continuously SOC mineralization in the incubation flasks. We calculated the delta SOC mineralized due to the change in soil moisture to assess the sensitivity of SOC mineralization to soil moisture regime. The second objective was to assess the sensitivity of SOC mineralization to temperature increase in carbon-storing soils. To do this, we also incubated replicates of the same soils in the laboratory at three temperatures (20, 28 and 35°C) in optimal moisture (pF 2.5). We calculated Q10 using the equal carbon respired method [START_REF] Conant | Sensitivity of organic matter decomposition to warming varies with its quality: temperature sensitivity of organic matter decomposition[END_REF] to assess the sensitivity of SOC mineralization to temperature in carbon-storing soils and baseline soils.

Materials and methods

Field site and soil sampling

This study focuses on two French long-term experiments (LTEs) where agroecological practices including conservation agriculture, organic agriculture and the application of organic waste products with contrasting biochemical quality have resulted in increased soil organic carbon contents and stocks relative to the baseline practices, i.e., conventional agriculture (Table . 1). This increase was named "additional soil organic carbon" (e.g., Bamière et al., 2022), and we used the term "carbon-storing soils" to refer to these soils that had stored additional carbon under agroecological practices.

La Cage LTE is conducted in Versailles (48°48'N,2°08'E). During the studied period , the mean annual temperature and precipitation were 11.6 °C and 633 mm respectively (Fig. 1).

The soil is a well-drained deep Luvisol (WRB, 2015). The experimental field is arranged in a randomized complete block design, divided into two blocks, themselves divided into four plots for each cropping system, and then into two subplots of 0.56 ha, so that wheat is present every year in one of the subplots [START_REF] Autret | Similar specific mineralization rates of organic carbon and nitrogen in incubated soils under contrasted arable cropping systems[END_REF]. A detailed presentation of crop rotations, soil management and fertilization were given by [START_REF] Autret | Alternative arable cropping systems: A key to increase soil organic carbon storage? Results from a 16 year field experiment[END_REF]. The 4 year's crop rotation mainly consisted of rapeseed (Brassica napus L.), winter wheat (Triticum aestivum L.), spring pea (Pisum sativum L.) and winter wheat. It differed in conservation agriculture (CA) and organic agriculture (ORG) for some years, with the replacement of rapeseed by maize (Zea mays L.) in CA or the introduction of alfalfa (Medicago sativa) in CA and ORG.

-CON-LC is characterized by a soil and crop management representative of the Paris Basin cereal production, with annual soil ploughing, the absence of organic amendment, a mineral N fertilization (average rate = 143 kg N ha -1 yr -1 ) and a systematic use of pesticides.

-CA includes a permanent soil cover, initially fescue (Festuca rubra) and since 2008 alfalfa, grown under the main crops, except pea. The soil is not tilled.

-ORG is characterized by an alfalfa-alfalfa-wheat-wheat rotation. No pesticides nor mineral fertilizers are used.

The QualiAgro LTE is located at Feucherolles, 20 km west of Versailles (48°52′N, 1°57′E) (Fig. 1). The soil is a Luvisol (WRB, 2015), cultivated for 21 years with a conventional wheat-maize rotation [START_REF] Peltre | RothC simulation of carbon accumulation in soil after repeated application of widely different organic amendments[END_REF]. The mean annual temperature and precipitation for the last 20 years are 11 °C and Research and Innovation since 1998, on which composts of organic waste products (OWPs) are applied every 2 years before tillage, at a dose equivalent to ~4 t C. ha -1 from 1998 to 2013 and ~2 t C. ha -1 from 2015 to 2020. The completely randomized block comprised an area of 6 ha with 40 plots of 450 m 2 each and 4 replicates per treatment. Since 2015, wheat and maize residues are buried in the soil. Three organic amendments are considered in this study and compared to a conventional agriculture treatment without organic inputs (CON-QA):

-Biowaste compost (BIOW): composting of the fermentable fraction of selectively collected household waste, mixed with green waste; -Municipal solid waste compost (MSW): composting of the residual fraction of household waste after selective collection of packaging; -Farmyard manure (FYM) which represents the reference amendment in the region. Soils from both LTEs were sampled at 30 ± 1 cm depth and stored in a chamber at 4°C after sieving to <4 mm a composite sample per plot.

Laboratory incubations

Surface soils (0-30cm) sieved to <4 mm was used for a microcosms incubation experiment to determine the effects of temperature (T°C), soil moisture (pF), and dry-wet-cycles (DWC) on SOC mineralization.

Before incubating the soil samples, a water content and matric potential response curve was established using Richard's presses. The water contents at saturation and corresponding to pF 1.5, pF 2.5 and pF 4.2 were determined for soils under each agricultural practice considered in this study (results presented in Table S1).

Soil cylinder construction and pre-incubation

Polyvinyl chloride (PVC) cylinders 5.7 cm in diameter and 4 cm in height with 2 mm perforations were used. A 50 μm mesh cloth at the bottom of the cylinder provided support for the soil while promoting gas exchange. Each cylinder was weighed empty and then with fresh soil equivalent to 100g dry soil.

Samples were then brought to a bulk density of 1.3 g.cm -3 with a hand press and mold. Knowing the initial water content, the samples were gradually brought to pF 2.5 adding water with a pasteur pipette.

Then the microcosms were mounted in 1L jars. The soil cylinders were placed on PVC supports and 15 mL of water was added to the bottom of the jars to stabilize the moisture. The jars were closed and the whole set was put in the incubator at 20°C for a one-week pre-incubation. Four replicates per agricultural practice and per moisture regime or temperature were prepared.

Experiment with soil moisture regime

To evaluate the effect of soil moisture regime on SOC mineralization, 4 soil moisture scenarios were performed and applied at 20°C: a continuously wet scenario at pF 1.5 (WET); a continuously moderated wet scenario at pF 2.5 (MWET); a continuously dry scenario at pF 4.2 (DRY) and a five subsequent dry-wet cycles scenario (DWC). To obtain the continuously moisture scenario after pre-incubation at pF 2.5, we added deionized water using a pasteur pipette to bring the replicates dedicated for WET scenario to pF 1.5. For replicates dedicated to MWET scenario (pF 2.5), we checked if the soil moisture corresponded to pF2.5. For the replicates dedicated to the DRY scenario, we introduced 50 g of silicagel during the pre-incubation phase to dry the soil. As soon as the moisture level corresponding to pF4.2 was reached, we removed the silica-gel from the jars. For DWC scenario, five dry-wet cycles were implemented during the experimental period (97 d). Each cycle contained two phases, 10 days of drying phase, followed by 7 days of moist phase after rewetting. The drying phase consisted of adding silica gel to the incubation jars to absorb the soil water and allow for a gradual soil drying. For this, perforated PVC cylinders containing 50 g of silica gel were added to each jar and changed after 6 days when the jars were vented after measuring the CO2 concentration in the headspace by micro-chromatography. A preliminary experiment had allowed to establish the amount of silica gel needed to dry the soil samples to pF 4.2 in 10 days. At the end of the drying period, rapid rewetting was performed by adding deionized water (the amount of water needed to reach pF 1.5) with a pasteur pipette. The same procedure was repeated for 5 times to simulate the dry-wet cycles in the soils. The four scenarios are schematized in Fig. 2.

Experiment with increased temperatures

To evaluate the effect of temperature on carbon mineralization, 3 temperature scenarios were compared: 20°C; 28°C and 35°C, using four replicates per scenario. The soil cylinders were kept at constant moisture (pF2.5). After the one-week pre-incubation the jars were placed either in an incubator set at 20°C, or in a thermostat chamber at 28°C and at 35°C and incubated for 92 days. We maintained constant soil moisture by weighing each sample after each CO2 measurement and adjusting the moisture content to the target mass.

Mineralization monitoring

Soil organic carbon (SOC) mineralization was monitored regularly for 97 days, with measurements on days 1, 3, 7, 13, 17, 24, 28, 34, 41, 45, 51, 59, 63, 69, 76, 80, 86 and 97 for each soil moisture scenario, and on days 1, 3, 7, 15, 22, 30, 37, 52, 65, 79 and 92 for each temperature scenario. SOC mineralization was measured non-destructively using a gas micro-chromatograph (μGC 490; Agilent Technologie; USA). The absolute amount of CO2 emitted is measured in parts per million (ppm). It is then converted to μg C-CO2 g -1 dry soil with the following formula:

(1), [START_REF] Védère | Spatial and temporal evolution of detritusphere hotspots at different soil moistures[END_REF], with CO2 (ppm): amount of CO2 emitted measured by gas micro-chromatograph; Mc: molar mass of carbon in g.mol -1 ; Vb: volume of the jar in L; VM: molar volume of gas in L.mol -1 and Msoil: mass of incubated soil in g.

Then the absolute amount of carbon mineralized was expressed per unit of SOC to obtain the specific SOC mineralization in μg C-CO2 /100 µg SOC, i.e., % SOC mineralized.

Soil moisture effect and sensitivity of SOC mineralization

The soil moisture effect on SOC mineralization was determined by calculating the delta SOC mineralized i.e., the difference in specific SOC mineralized between the high moisture (WET), lowmoisture (DRY) and dry-wet-cycle (DWC) with specific SOC mineralized under optimal moisture (MWET). A negative value indicates an inhibition of SOC mineralization, while a positive value indicates a stimulation of SOC mineralization. Furthermore, to assess the sensitivity of carbon-storing soils to the soil moisture regime, we compared the delta SOC mineralized under WET, DRY, and DWC between baseline soils (CON-LC and CON-QA for La Cage and QualiAgro LTEs, respectively) and soils under agroecological practices, i.e., CA, ORG, MSW, FYM and BIOW.

Temperature sensitivity of SOC mineralization

The Q10 was used in this study to estimate the sensitivity of organic carbon mineralization under different agricultural practices to temperature increase. The method for calculating Q10 at equal carbon mineralized (Q10-q) described by [START_REF] Conant | Sensitivity of organic matter decomposition to warming varies with its quality: temperature sensitivity of organic matter decomposition[END_REF] was carried out in this study. This method involves determining the amount of time needed for a given amount of SOC to be respired at a given temperature. The time required to respire a given amount of SOC at two temperatures is then used to calculate a Q10 value. This method of analysis eliminates the problem of fixed incubation duration leading to comparison of different SOC pools, which confounds characterization of temperature sensitivity [START_REF] Reichstein | Temperature dependence of carbon mineralisation: conclusions from a long-term incubation of subalpine soil samples[END_REF], 2005[START_REF] Leifeld | The Temperature Response of CO2 Production from Bulk Soils and Soil Fractions is Related to Soil Organic Matter Quality[END_REF]. Therefore, the Q10 was calculated for 1% and 3% of respired SOC. The Q10 is calculated for 2 temperature ranges (20-28°C and 28-35°C) with the following formula:

(2) [START_REF] Conant | Sensitivity of organic matter decomposition to warming varies with its quality: temperature sensitivity of organic matter decomposition[END_REF], with t1 and t2 the time required to respired the same amount of carbon at low temperature (T1) and high temperature (T2).

Statistical analysis

All data were tested for normality and homogeneity of variance. Log-transformation was applied, if the transformation improved the normality and variance substantially. A one-way ANOVA with Tukey's test was used to detect the differences in Q10, the delta SOC mineralized and the amount of SOC mineralized among soil moisture and temperature scenarios within each field site, then across all scenarios in both sites. When the normality or homogeneity of the data was not confirmed, we applied the non-parametric Kruskal-Wallis test. All statistical analyses were completed in R (version 4.0.2).

Results

Soil carbon mineralization under diverse moisture regimes, temperatures and agricultural practices

Increasing soil moisture significantly affected SOC mineralization in the La Cage and QualiAgro LTEs (Table. 2; p<0.05). Also, SOC mineralization followed a polynomial function with the soil matric potential (Fig. 3). Carbon mineralization under soil moisture regime was in this order: WET=MWET>DWC>DRY. To better observe the effect of soil moisture regime on carbon mineralization, we considered the MWET as the moisture control in order to calculate the delta SOC mineralization induced by the change in soil moisture. The delta SOC mineralized in the WET scenario compared to MWET was on average +0.07 to +1.10 % of SOC (Fig. 6), showing a lack of or a weak stimulation of SOC mineralization with increasing soil moisture above pF 2.5 (Fig. 6). The delta SOC mineralized was significantly larger at La Cage (CON-LC: 1.10 ± 0.45; ORG: 0.89 ± 0.76 and CA: 0.69 ± 0.52% of SOC) than at QualiAgro (CON-QA: 0.12 ± 0.43; MSW: 0.24 ± 0.24; FYM: 0.07 ± 0.12 and BIOW: 0.30 ± 0.10% of SOC) (p<0.05; Fig. 6). Increasing the soil moisture from pF2.5 to pF1.5 did not significantly increase carbon mineralization at QualiAgro. In contrast, the DWC and DRY scenarios inhibited carbon mineralization relative to MWET scenario (Fig. 6). SOC mineralization in DWC decreased on average between 1.15 to 2.16 % of SOC and in DRY scenario on average between 2.35 to 2.97 % SOC (Fig. 6). This result suggests that multiple dry-wet cycles (DWC) did not stimulate carbon mineralization relative to the optimum soil moisture (MWET), whereas multiple DWC stimulated carbon mineralization relative to the low moisture scenario (DRY).

Increasing temperature also significantly affected cumulative SOC mineralization in both LTEs (Table 3; p<0.05). Except the CON-LC where no significant difference was found between mineralization at 20 and 28°C, the soils under other agricultural practices showed significant differences between 20, 28 and 35°C (Table . 3; p<0.05).

At QualiAgro, we observed a significantly larger specific SOC mineralization in the reference plots (CON-QA) compared to the plots that received organic waste products (CON-QA>BIOW=FYM =MSW), whatever the soil moisture regime and temperature (Fig. 5, Fig. 7c, Fig. 7d), indicating less stabilization of SOC in CON-QA. However, at La Cage, regardless of soil moisture regime and temperature (Fig. 4, Fig. 7a, Fig. 7b), there was no difference in specific mineralization rates between the treatments (ORG, CA) relative to the conventional agriculture (CON-LC).

In addition, the effect of agricultural practices in both LTEs reported in Table 2, indicates that CON-LC baseline practice at La Cage and the CON-QA baseline practice at QualiAgro, not only had similar SOC contents (9.82 ± 0.48 g C.kg -1 and 10.11 ± 1.12 g C.kg -1 respectively) but also had the same percentage of mineralized SOC whatever the soil moisture regime. Moreover, the percentage of SOC mineralized under agroecological practices at La Cage (CA: 2.24 to 5.61 % SOC and ORG: 2.56 to 6.38 % SOC) was significantly larger than under organic waste products application (MSW: 1.46 to 4.09 % SOC; FYM: 1.48 to 4.06 % SOC and BIOW: 1.23 to 3.89 % SOC) (Table . 2). Also, when comparing the agroecological practices in both LTEs, we found that the agroecological practices at La Cage (ORG and CA) generated more SOC loss than agroecological practices at QualiAgro (BIOW, MSW, and FYM) after 3 months of incubation at diverse temperatures (Table . 3).

Sensitivity of SOC mineralization to soil moisture regime

The analysis of variance (ANOVA) revealed that within each LTE, we observed no significant difference between the carbon-storing soils (CA, ORG, MSW, FYM and BIOW) and their respective baseline soils (CON-LC and CON-QA) in the delta SOC mineralized whatever the soil moisture regime (Fig. 6). These results indicate that the carbon-storing soils had a similar sensitivity to the soil moisture regime as the baseline ones.

Sensitivity of SOC mineralization to temperature increase

The mean Q10 corresponding to equal 1% mineralized SOC at QualiAgro ranged between 1.90 to 2.53 for the 20-28°C temperature range, compared to 1.38 to 1.55 at La Cage (Fig. 8). These Q10 did not change significantly from 1% to 3% mineralized SOC (Table S2). Furthermore, the Q10 for the temperature range 20-28°C was larger for the BIOW plot than for the CON-LC plot (Fig. 8; p<0.05), but no significant difference was found between agroecological practices (CA, ORG, MSW, FYM and FYM). For the temperature range 28-35°C we did not observe any significant difference between agricultural practices (Fig. 8, p>0.05). The SOC under OWPs practices (BIOW, FYM and MSW) had the same sensitivity of SOC mineralization to temperature as their baseline practice (CON-QA), as did the CA and ORG with their baseline practice (CON-LC). These results indicate that carbon-storing soils had a similar sensitivity of SOC mineralization to temperature increase as the baseline soils.

Discussion

Soil carbon mineralization under soil moisture regime and temperature

The interaction of soil microbes with their physical environment affects their ability to mineralize SOC into CO2. One of these environmental factors is soil moisture [START_REF] Cook | Relationships between soil respiration and soil moisture[END_REF].

The results of the present study indicated that SOC mineralization followed a polynomial function with soil matric potential (Fig. 3). [START_REF] Manzoni | Responses of soil microbial communities to water stress: results from a meta-analysis[END_REF] and [START_REF] Moyano | The moisture response of soil heterotrophic respiration: interaction with soil properties[END_REF] in their meta-analyses found for their dataset that the soil heterotrophic respiration (mineralization) decreased steadily when the matric potential became more negative, i.e., for high pF values. This result agrees with ours, where SOC mineralization decreased as we reach higher pFs (pF1.5 ³ pF2.5>pF4.2), i.e., lower matric potentials (-0.015 MPa ³ -0.033 MPa > -1.6 MPa). Such a decrease in respiration is associated with a reduction in solute diffusivity in mineral soils at low soil moisture (pF4.2), because water-filled capillaries become disconnected as the soil becomes drier [START_REF] Moldrup | Tortuosity, Diffusivity, and Permeability in the Soil Liquid and Gaseous Phases[END_REF]. Therefore, solute diffusivity (and thus substrate and nutrient bioavailability) may be the most limiting factor under dry conditions [START_REF] Skopp | Steady-State Aerobic Microbial Activity as a Function of Soil Water Content[END_REF][START_REF] Schjønning | Linking Soil Microbial Activity to Water-and Air-Phase Contents and Diffusivities[END_REF][START_REF] Or | Physical constraints affecting bacterial habitats and activity in unsaturated porous media -a review[END_REF], Moyano et al., 2013). Conversely, increasing soil moisture increases solute diffusion rates in soils, enhancing the accessibility of substrates by soil microorganisms [START_REF] Monard | Impact of soil matric potential on the fine-scale spatial distribution and activity of specific microbial degrader communities[END_REF], Moyano et al., 2013[START_REF] Manzoni | Invariant soil water potential at zero microbial respiration explained by hydrological discontinuity in dry soils[END_REF]. Our results are consistent with the hypothesis that higher soil moisture increase the accessibility of SOC to decomposers and thus lead to more SOC mineralization than low moisture contents, provided the soil moisture remains well below saturation, where oxygen diffusion rates becomes a limiting factor for heterotrophic respiration (Cook andOrchard, 2008, Moyano et al., 2013).

In addition, the percentage of SOC mineralized at pF1.5 (WET) after 97 days are statistically identical to the percentage at pF2.5 (MWET) (Table . 2), implying that the increase in soil moisture beyond pF2.5 did not result in significant further CO2 emission to the atmosphere. This result suggests that at pF2.5, the mineralization optimum was reached, and increasing soil moisture to pF1.5 does not yet alter SOC mineralization due to slight anoxia. Because the diffusion rate of oxygen through water is much lower than through air [START_REF] Cook | Oxygen Transport to Plant Roots: Modeling for Physical Understanding of Soil Aeration[END_REF], the metabolic activity of aerobic organisms also decreases as soil pore space fills with water and approaches saturation levels [START_REF] Franzluebbers | Microbial activity in response to water-filled pore space of variably eroded southern Piedmont soils[END_REF]. [START_REF] Curtin | Temperature and Moisture Effects on Microbial Biomass and Soil Organic Matter Mineralization[END_REF] showed that the optimum matric potential for mineralization can be as high as -0.005 MPa; while the meta-analysis by [START_REF] Moyano | The moisture response of soil heterotrophic respiration: interaction with soil properties[END_REF] found an optimum matric potential for mineralization of -0.001 MPa corresponding to pF 1. In our case, we observed no decrease in SOC mineralization in the wettest state, i.e., -0.015 MPa (pF1.5). One possible explanation is that the soil was sieved and repacked and the cores had a high inter-aggregate macroporosity favorable to air circulation, preventing anoxia, or that at pF1.5 a slight anoxia at compensated for a better substrate diffusion compared to pF2.5 resulting in similar SOC mineralization rates at pF1.5 and 2.5.

The dry-wet cycles scenario tested in this study did not induce substantial delta SOC mineralization when compared to optimum moisture (MWET), but did if compared to low moisture (DRY) scenario (Table . 2). The drying phase induced low carbon mineralization rates as a result of reduced microbial activity (Franzluebbers et al., 1994, Pulleman and[START_REF] Pulleman | Microbial C and N transformations during drying and rewetting of coniferous forest floor material[END_REF] and decreased accessibility of the organic substrate to decomposers [START_REF] Manzoni | Responses of soil microbial communities to water stress: results from a meta-analysis[END_REF]. This low mineralization was not compensated by the flush observed after rewetting [START_REF] Mikha | Carbon and nitrogen mineralization as affected by drying and wetting cycles[END_REF][START_REF] Yemadje | Wetting-drying cycles do not increase organic carbon and nitrogen mineralization in soils with straw amendment[END_REF]. Zhang et al. (2020) in their meta-analysis showed that, on average, dry-wet-cycles stimulated soil carbon mineralization by 72% relative to mineralization at low-moisture and inhibited carbon mineralization by 25% relative to incubation at high-moisture which is consistent with our results. Relevant controls are indeed needed when evaluating the importance of the Birch effect (Kpemoua et al., 2023).

Temperature generally had a positive effect on cumulative SOC mineralization, which was manifested by larger CO2 emissions at higher temperatures (Table . 3). Similar increases in SOC mineralization with increasing temperature have previously been observed in many other studies (e.g., [START_REF] Conant | Sensitivity of organic matter decomposition to warming varies with its quality: temperature sensitivity of organic matter decomposition[END_REF][START_REF] Ghosh | Temperature sensitivity of soil organic carbon decomposition as affected by long-term fertilization under a soybean based cropping system in a sub-tropical Alfisol[END_REF][START_REF] Ghimire | Long-term Management Effects and Temperature Sensitivity of Soil Organic Carbon in Grassland and Agricultural Soils[END_REF][START_REF] Fu | Increasing temperature can modify the effect of straw mulching on soil C fractions, soil respiration, and microbial community composition[END_REF]. Temperature rise affects most of the physiological activity of microbial cell and accelerates residue decomposition [START_REF] Reichstein | Temperature dependence of carbon mineralisation: conclusions from a long-term incubation of subalpine soil samples[END_REF][START_REF] Yun | Effects of Temperature and Moisture on Soil Organic Carbon Mineralization[END_REF].

Soil carbon mineralization under agricultural practices

We observed no significant difference in specific SOC mineralization rate between conventional agriculture (CON-LC) and agroecological practices (ORG and CA) at La Cage, whatever soil moisture regime and whatever the temperature. This result supports previous results obtained by [START_REF] Autret | Similar specific mineralization rates of organic carbon and nitrogen in incubated soils under contrasted arable cropping systems[END_REF] on the same LTE, where soil incubation at pF 3 reported similar specific carbon mineralization between alternative practices (conservation, organic and low input agriculture) and conventional agriculture. Some studies (Dendooven et al., 2012, Ussiri and[START_REF] Ussiri | Long-term tillage effects on soil carbon storage and carbon dioxide emissions in continuous corn cropping system from an alfisol in Ohio[END_REF] did measure lower SOC mineralization rates under no tillage, ascribed to less favorable local climatic conditions and better physical protection of organic matter in the soil structure [START_REF] Balesdent | Relationship of soil organic matter dynamics to physical protection and tillage[END_REF]. Similarly to our results, Zhang et al. (2020b), showed that tillage had no influence on specific C mineralization, while the interaction of tillage and depth had a significant effect. In the present study, we incubated samples from the surface to a depth equivalent to the soil mass of 4300 kg. ha -1 (corresponding to 30 ± 1 cm), regardless of agricultural practice. No-tillage typically results in carbon stratification in the soil (e.g., [START_REF] Zhao | Stratification and Storage of Soil Organic Carbon and Nitrogen as Affected by Tillage Practices in the North China Plain[END_REF]. We hypothesize that the physical protection of carbon is increased in the first 10 cm under no-tillage, but that considering a 0-30 cm soil layer cancels out this effect. Another hypothesis is that sieving the soils to <4 mm before incubation would have removed the physical protection occurring in macroaggregates in the CA soils.

At the QualiAgro LTE, the application of organic waste products (MSW, FYM and BIOW) reduced significantly SOC mineralization, compared to conventional agriculture without organic inputs (CON-QA) whatever soil moisture regime and whatever the temperature. These results can be explained by the fact that composts and farmyard manure have previously undergone a period of decomposition and therefore mineralize less rapidly than crop residues that contain higher proportions of easily mineralizable organic carbon [START_REF] Chodak | Effect of temperature on the mineralization of C and N of fresh and mature compost in sandy material[END_REF][START_REF] Leifeld | Biological activity and organic matter mineralization of soils amended with biowaste composts[END_REF]. The physico-chemical composition of the organic waste products (OWPs) has been analyzed by [START_REF] Peltre | RothC simulation of carbon accumulation in soil after repeated application of widely different organic amendments[END_REF] at the QualiAgro LTE. They show that residual municipal solid waste compost (MSW) is the OWP with the least lignin, and biowaste compost (BIOW) the one with the most. Lignin is known to be a compound resistant to microbial degradation at a year to decadal timescale, which would explain why BIOW tends to mineralize less SOC than MSW and FYM (Table . 2). In the same way, an incubation study of substrates of different quality showed that manure is less biodegradable than crop residues such as straw and suggests that this is due to a more stable chemical quality, as manure is already biodegraded during animal digestion and storage [START_REF] Benbi | Effects of Temperature, Moisture, and Chemical Composition of Organic Substrates on C Mineralization in Soils[END_REF]. [START_REF] Fortuna | Optimizing nutrient availability and potential carbon sequestration in an agroecosystem[END_REF] found that easily mineralizable C constituted less of the SOC in soils fertilized with compost as compared with soils that were fertilized with mineral N. Other studies have also shown that the application of organic amendments that include composts and manure have the potential to reduce CO2 emission per unit of SOC since these best management practices contribute to carbon stabilization [START_REF] Dou | Sensitivity of Labile Soil Organic Carbon to Tillage in Wheat-Based Cropping Systems[END_REF][START_REF] Bhowmik | Potential carbon sequestration and nitrogen cycling in long-term organic management systems[END_REF].

Overall, the cumulative specific SOC mineralization in our incubations was much less for soils that received organic waste products plus crop residues (MSW, FYM, BIOW) than for soils receiving only crop residues (CON-QA, CON, ORG, CA) (Table . 2 andTable. 3). This result thus confirms that soils that received organic waste were depleted in labile SOC compared to soils that received only crop residues. In contrast to our results, [START_REF] Obriot | Epandage de produits résiduaires organiques et fonctionnement biologique des sols: De la quantification des impacts sur les cycles carbone et azote à l'evaluation multicritère de la pratique à l'échelle de la parcelle[END_REF] had observed similar proportions of soil organic carbon mineralization between the conventional system without organic inputs and plots that had received 7 successive applications of organic waste products. We assume that the long-term application (11 successive applications in our study) of organic waste products led to a greater stabilization of carbon in the soil, which makes it less biodegraded by microorganisms.

Sensitivity of SOC mineralization to soil moisture regime

The results of the delta SOC mineralized due to soil moisture regime change indicated no significant difference between carbon-storing soils (CA and ORG) and the baseline soil (CON-LC) at La Cage (Fig. 6). These results suggest that the decomposability of organic carbon is similar under the contrasted agricultural practices whatever the soil moisture regime, or the environmental conditions for decomposers are the same [START_REF] Autret | Similar specific mineralization rates of organic carbon and nitrogen in incubated soils under contrasted arable cropping systems[END_REF]. Regarding the decomposability of the organic carbon, it is likely that the nature of organic inputs to soil is rather similar as soil organic matter exhibits similar C/N ratio (Table . 1) and wheat is cropped one year out of two in all plots [START_REF] Autret | Alternative arable cropping systems: A key to increase soil organic carbon storage? Results from a 16 year field experiment[END_REF]. The sensitivity of SOC mineralization to soil moisture regime is not affected by the cropping system nor tillage practice in this LTE.

At the QualiAgro LTE, the delta SOC mineralized due to changes in soil moisture regime also showed no significant difference between the baseline soil (CON-QA) and carbon-storing soils (MSW, BIOW and FYM). This result indicates that although the SOC in the baseline soil is more biodegradable on average, it has similar sensitivity of carbon mineralization to soil moisture regime than the SOC of carbon storing soils. The results in both experiments are contrary to our initial hypothesis concerning the higher sensitivity of carbon mineralization to soil moisture regime in carbon-storing soils relative to their baseline. Therefore, in a drier, wetter, or dry-wet conditions, the agroecological practices implemented in both LTEs are still beneficial for climate change mitigation.

Sensitivity of SOC mineralization to temperature increase

We found that the Q10 values (20-28°C: 1.38, 1.52 and 1.55 respectively for CON-LC, CA and ORG) were not significantly different under the contrasted agricultural practices at La Cage (Fig. 8).

These results are contrary to those reported by [START_REF] Parihar | Dependence of temperature sensitivity of soil organic carbon decomposition on nutrient management options under conservation agriculture in a sub-tropical Inceptisol[END_REF], who found the extent of increase in mineralization with temperature elevation to be higher under conventional agriculture (high Q10) than under conservation agriculture (low Q10) whatever the soil depth. They justified their results by the high physical protection of the SOC under conservation agriculture limiting its accessibility to soil microbial degradation [START_REF] Six | Soil organic matter, biota and aggregation in temperateand tropical soils -Effects of no-tillage[END_REF]. However, [START_REF] Autret | Alternative arable cropping systems: A key to increase soil organic carbon storage? Results from a 16 year field experiment[END_REF][START_REF] Autret | Similar specific mineralization rates of organic carbon and nitrogen in incubated soils under contrasted arable cropping systems[END_REF] had found in the La Cage LTE that the increase in carbon stocks in conservation agriculture was mainly due to the carbon input rather than to the expected physical protection. This result confirms ours, leading us to assume that neither tillage practice nor cropping systems affected the sensitivity of SOC mineralization to temperature increase. So, the carbon-storing soils (CA, ORG) at La Cage had a similar sensitivity of SOC mineralization to temperature increases as the baseline soil (CON-LC).

At QualiAgro, the Q10 value (Fig. 8) also did not reveal any significant difference between the agricultural practices (20-28°C: 1.90; 2.44; 2.37 and 2.52 for CON-QA, MSW, FYM and BIOW respectively). These results suggest that the carbon-storing soils (MSW, FYM, BIOW) at QualiAgro and the baseline soil (CON-QA) had the same sensitivity of SOC mineralization to temperature increase. This result contrasts with [START_REF] Benbi | Effects of Temperature, Moisture, and Chemical Composition of Organic Substrates on C Mineralization in Soils[END_REF] who showed that temperature sensitivity depends on the quality of the substrate, and that manure is more sensitive (Q10 = 3) than green waste which is easily decomposable (Q10 = 2.5). However, they measured the mineralization of manure and green waste directly and not the soils that were amended with manure or green waste as we did.

In summary, the Q10 values obtained in the present study were generally close to the values reported in the meta-analysis by [START_REF] Hamdi | Synthesis analysis of the temperature sensitivity of soil respiration from laboratory studies in relation to incubation methods and soil conditions[END_REF] for agricultural soils. No differences were observed between the agricultural practice in both LTEs on their temperature sensitivity, although slightly higher Q10 values were observed at QualiAgro, especially in OWPs practices than in conservation and organic agriculture. To test the robustness of this result, we also calculated the Q10 calculated for 3% of SOC being mineralized in the incubations. We found it was not significantly different, neither across agricultural practices nor from the Q10 calculated for 1% SOC mineralized (Table S2). This result suggests that the carbon pool mineralized during 3 months incubation had similar quality and sensitivity whatever the agricultural practice.

Conclusion

While a diversity of agroecological management options allow to store additional organic carbon in soil, a remaining question concerns the permanence and vulnerability of this stored carbon. This study examined the response of soil organic carbon stored by agroecological practices to climate change. We found that the heterotrophic respiration of "carbon-storing soils" had similar sensitivities compared to their baseline counterparts, regarding soil moisture regime changes and temperature increases. Hence, the implementation of these agroecological practices appears beneficial for climate change mitigation, even in the context of extreme climatic events.

Here were assessed the sensitivity of soils managed under different agricultural practices to climate change, i.e., according to IPCC (2001) the degree to which a system or species is affected, either adversely or beneficially, by climate variability or change. IPCC considers that the vulnerability to climate change is wider, as it encompasses also the capacity to adapt to climate change. In further studies, other aspects could be accounted for, such as the soils capacity to retain water, which our results suggest to be increased in the carbon-storing soils, the changes moisture and temperature regime consequent to the presence of a mulch, such as in conservation agriculture, or to changes in albedo, consequent to OWPs applications.

List of tables

Table 1 Physical and chemical properties of the soil at "La Cage" (layer 0-25 cm) and "QualiAgro (layer 0-29 cm)" measured at the start of experiment in 1998 [START_REF] Peltre | RothC simulation of carbon accumulation in soil after repeated application of widely different organic amendments[END_REF][START_REF] Autret | Alternative arable cropping systems: A key to increase soil organic carbon storage? Results from a 16 year field experiment[END_REF] to 2020 at QualiAgro and La Cage (layer 0-30 cm). CON-LC (conventional agriculture); ORG (organic agriculture); CA (conservation agriculture); CON-QA (conventional agriculture without organic inputs); MSW (residual municipal solid waste compost); FYM (farmyard manure) and BIOW (biowaste compost).

Table 2 Soil organic carbon content according to agricultural practices, percentage of SOC mineralized according to soil moisture scenario: continuously wet scenario at pF1.5 (WET); continuously moderated wet scenario at pF2.5 (MWET); continuously dry scenario at pF4.2 (DRY) and Dry-wet cycles scenario (DWC). SOC: soil organic carbon; CON-LC: conventional agriculture; ORG: organic agriculture; CA: conservation agriculture; CON-QA: conventional agriculture without organic waste products; MSW: residual municipal solid waste compost; FYM: farmyard manure; BIOW: biowaste compost. NS: nonsignificant; ANOVA: analysis of the variance. (¡): p-value of soil moisture regime effect on SOC mineralization and (^): p-value of agricultural practices effect. Uppercase letters represent differences between soil moisture scenarios and lowercase letters represent differences between agricultural practices. A significant difference is obtained for p<0.05.

Table 3 Percentage of SOC mineralized according to temperature increasing. (¡): p-value of incubation temperature effect on SOC mineralization and (^): p-value of agricultural practices effect. Uppercase letters represent differences between temperature scenarios and lowercase letters represent differences between agricultural practices. A significant difference is obtained for p<0.05. Table . 1 Physical and chemical properties of the soil at "La Cage" (layer 0-25 cm) and "QualiAgro (layer 0-29 cm)" measured at the start of experiment in 1998 [START_REF] Peltre | RothC simulation of carbon accumulation in soil after repeated application of widely different organic amendments[END_REF][START_REF] Autret | Alternative arable cropping systems: A key to increase soil organic carbon storage? Results from a 16 year field experiment[END_REF] to 2020 at QualiAgro and La Cage (layer 0-30 cm). CON-LC (conventional agriculture); ORG (organic agriculture); CA (conservation agriculture); CON-QA (conventional agriculture without organic inputs); MSW (residual municipal solid waste compost); FYM (farmyard manure) and BIOW (biowaste compost).

List of figures

Table. 2

Soil organic carbon content according to agricultural practices, percentage of SOC mineralized according to soil moisture scenario: continuously wet scenario at pF1.5 (WET); continuously moderated wet scenario at pF2.5 (MWET); continuously dry scenario at pF4.2 (DRY) and Dry-wet cycles scenario (DWC). SOC: soil organic carbon; CON-LC: conventional agriculture; ORG: organic agriculture; CA: conservation agriculture; CON-QA: conventional agriculture without organic waste products; MSW: residual municipal solid waste compost; FYM: farmyard manure; BIOW: biowaste compost. NS: nonsignificant; ANOVA: analysis of the variance. (¡): p-value of soil moisture regime effect on SOC mineralization and (^): p-value of agricultural practices effect. Uppercase letters represent differences between soil moisture scenarios and lowercase letters represent differences between agricultural practices. A significant difference is obtained for p<0.05.

Fig. 1 :

 1 Fig. 1: Map of the Yvelines department (France) showing the location of the two field sites.

Fig. 2 .

 2 Fig. 2. A schematic diagram of the experimental design showing constant moisture scenarios (WET,

Fig. 3 .

 3 Fig. 3. Relationship between cumulative carbon mineralization at the end of incubation (97 days) and

Fig. 4 .

 4 Fig. 4. Continuously soil moisture (WET, DRY, MWET) and dry-wet-cycles (DWC) effect on specific

Fig. 5 .

 5 Fig. 5. Continuously soil moisture (WET, DRY, MWET) and dry-wet-cycles (DWC) effect on specific

Fig. 6 .

 6 Fig. 6. Delta SOC mineralized under WET, DWC and DRY scenarios relative to MWET scenario. Error

Fig. 7 .

 7 Fig. 7. Temperature effect on specific cumulative soil organic carbon (SOC) mineralization under 28°C

Fig. 8 .

 8 Fig. 8. Q10 values distribution according to agricultural practices. On the left the Q10 for the temperature

Fig. 1 :

 1 Fig. 1: Map of the Yvelines department (France) showing the location of the two field sites.

Fig. 2 .

 2 Fig. 2. A schematic diagram of the experimental design showing constant moisture scenarios (WET,

Fig. 3 .

 3 Fig. 3. Relationship between cumulative carbon mineralization at the end of incubation (97 days) and

Fig. 4 .

 4 Fig. 4. Continuously soil moisture (WET, DRY, MWET) and dry-wet-cycles (DWC) effect on specific

Fig. 5 .

 5 Fig. 5. Continuously soil moisture (WET, DRY, MWET) and dry-wet-cycles (DWC) effect on specific

Fig. 6 .

 6 Fig. 6. Delta SOC mineralized under WET, DWC and DRY scenarios relative to MWET scenario. Error

Fig. 7 .

 7 Fig. 7. Temperature effect on specific cumulative soil organic carbon (SOC) mineralization under 28°C

Fig. 8 .

 8 Fig. 8. Q10 values distribution according to agricultural practices. On the left the Q10 for the temperature

Table . 3

 . Percentage of SOC mineralized according to temperature increasing. (¡): p-value of incubation temperature effect on SOC mineralization and (^): p-value of agricultural practices effect. Uppercase letters represent differences between temperature scenarios and lowercase letters represent differences between agricultural practices. A significant difference is obtained for p<0.05.

						SOC				% SOC mineralized
		LTE	Agricultural	g.kg -1	WET		MWET		DWC	DRY	p-value (¡) ANOVA
			practices
			CON-LC	9.82 ± 0.48 d	5.53 ± 0.44 abA	4.43 ± 0.74 aA	3.27± 0.47 abC	1.93 ± 0.26 abD
		La Cage		ORG CA	10.39 ± 0.42 d 13.30 ± 1.05 b	6.38 ± 0.52 aA 5.61 ± 0.88 abA	5.57 ± 0.31 aA 4.67 ± 0.26 aA	3.33 ± 0.43 aC 3.28 ± 0.67 abB	2.56 ± 0.33 aD 2.24 ± 0.26 aC	p<0.05
			CON-QA	10.11 ± 1.12 d	4.79 ± 0.26 bcA 4.92 ± 0.22 abA	3.22 ± 0.26 abB	1.95 ± 0.49 bcC
		QualiAgro	MSW	13.38 ± 0.45 b	4.09 ± 0.48 cA	3.86 ± 0.33 bA	2.55 ± 0.08 bB	1.46 ± 0.19 cC
				FYM	14.31 ± 0.63 ab	4.06 ± 0.24 cA	3.99 ± 0.24 bA	2.65 ± 0.11 bB	1.48 ± 0.17 cC	p<0.05
				BIOW	15.81 ± 1.23 a	3.89 ± 0.25 cA	3.58 ± 0.18 bA	2.43 ± 0.19 bB	1.23 ± 0.08 cC
		p-value (^)		p<0.05	p<0.05		p<0.05	p<0.05	p<0.05
	LTEs	Agricultural	Tillage ANOVA	Block	Clay	Silt	Sand	C	N	C/N	pH H20	Exchangeable cations (Cmol + .kg -1 )	CEC
		practices	modality				g.kg -1					K +	Na +	Ca 2+	Mg 2+	Cmol + .
																	kg -1
	La Cage 1998			-				9.85	0.99	9.95	7.37	-	-	-	-	11.52
	La Cage	CON-LC			1				10.05	0.88	11.55	7.20	0.546	0.039	13.80	1.78	11.10
	2020		Tillage	2	170	560	270	9.25	0.91	10.15	7.15	0.584	0.039	13.30	1.75	9.85
		ORG			1				10.60	0.99	10.65	6.90	0.395	0.039	13.40	1.97	10.10
					2				10.05	0.95	10.60	7.20	0.474	0.039	13.00	2.01	9.80
		CA	No		1				13.55	1.32	10.25	7.15	0.707	0.039	13.95	2.28	10.95
			Tillage	2				12.05	1.06	11.50	7.15	0.679	0.039	13.70	2.33	10.75
	QualiAgro 1998			-				10.39	1.12	9.25	6.80	0.442	0.021	8.298	0.494	9.71
	QualiAgro	CON-QA			-				10.11	1.04	9.65	6.74	0.270	0.032	7.67	0.455	8.25
	2020	MSW	Tillage	-	156	778	65	13.38	1.31	10.25	7.77	0.375	0.053	10.19	0.533	10.75
		FYM			-				14.30	1.41	10.16	7.31	0.922	0.039	8.65	0.950	10.27
		BIOW			-				16.62	1.66 % SOC mineralized 9.12 7.71	0.610	0.04	10.86	0.650	11.77
		LTEs	Agricultural practices		20°C			28°C		35°C		p-value (¡) ANOVA
				CON-LC		4.28 ± 0.72 aA		5.90 ± 1.07 bA		7.96 ± 0.79 bB
		La Cage			ORG CA		5.37 ± 0.31 aA 5.01 ± 0.24 aA		7.37 ± 0.52 abB 7.26 ± 0.33 abB		8.35 ± 0.36 bC 8.55 ± 0.52 bC		p<0.05
				CON-QA		4.77 ± 0.48 aA		7.89 ± 0.91 aB		9.54 ± 0.46 aC
		QualiAgro		MSW		3.74 ± 0.33 bA		6.01 ± 0.67 bB		7.91 ± 0.29 bC
				FYM		3.87 ± 0.23 bA		6.13 ± 0.42 bB		8.40 ± 0.44 bC		p<0.05
				BIOW		3.47 ± 0.18 bA		5.87 ± 0.51 bB		7.52 ± 0.41 bC
		p-value (^) ANOVA/Kruskal test		p=0.02			p<0.05		p<0.05

Acknowledgements: This research has been supported by the French Agence Nationale de la Recherche (StoreSoilC project, grant ANR-17-CE32-0005). C. Chenu also thanks the CLand programme (ANR-16-CONV-0003). T. P. I. Kpemoua acknowledges the support of ADEME. The QualiAgro field experiment forms part of the SOERE-PRO (network of long-term experiments dedicated to the study of impacts of organic waste product recycling) integrated as a service of the ''Investment for future'' infrastructure AnaEE-France, overseen by the French National Research Agency (ANR-11-INBS-0001)'. The QualiAgro site is conducted within a collaboration between INRAE and Veolia. The authors also thank the personnel from Unité expérimentale Versailles Grignon INRAE and Michel Bertrand for the maintenance of and access to the La Cage long term experiment.

Supplementary materials for:

Are carbon-storing soils more sensitive to climate change? A laboratory evaluation for agricultural temperate soils Tchodjowiè P. I. Kpemoua 1,3


Subhasish Basak

Sébastien Petit

Julien Bect

Emmanuel Vazquez

Numerical issues in maximum likelihood parameter estimation for Gaussian process interpolation

Keywords: Gaussian process, Maximum likelihood estimation, Optimization

This article investigates the origin of numerical issues in maximum likelihood parameter estimation for Gaussian process (GP) interpolation and investigates simple but effective strategies for improving commonly used opensource software implementations. This work targets a basic problem but a host of studies, particularly in the literature of Bayesian optimization, rely on off-theshelf GP implementations. For the conclusions of these studies to be reliable and reproducible, robust GP implementations are critical.

Introduction

Gaussian process (GP) regression and interpolation (see, e.g., [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], also known as kriging (see, e.g., [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF], has gained significant popularity in statistics and machine learning as a non-parametric Bayesian approach for the prediction of unknown functions. The need for function prediction arises not only in supervised learning tasks, but also for building fast surrogates of time-consuming computations, e.g., in the assessment of the performance of a learning algorithm as a function of tuning parameters or, more generally, in the design and analysis computer experiments [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. The interest for GPs has also risen considerably due to the development of Bayesian optimization [START_REF] Mockus | On Bayesian methods for seeking the extremum[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Emmerich | Single-and multiobjective evolutionary optimization assisted by Gaussian random field metamodels[END_REF][START_REF] Srinivas | Gaussian process optimization in the bandit setting: no regret and experimental design[END_REF].

This context has fostered the development of a fairly large number of open-source packages to facilitate the use of GPs. Some of the popular choices are the Python modules scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], GPy (Sheffield machine learning group, 2020), GPflow [START_REF] Matthews | GPflow: A Gaussian process library using TensorFlow[END_REF], GPyTorch [START_REF] Gardner | GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration[END_REF], OpenTURNS [START_REF] Baudin | OpenTURNS: an industrial software for uncertainty quantification in simulation[END_REF]; the R package DiceKriging [START_REF] Roustant | DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF]; and the Matlab/GNU Octave toolboxes GPML [START_REF] Rasmussen | Gaussian processes for machine learning (GPML) toolbox[END_REF], STK (Bect et al., 2021) and GPstuff [START_REF] Vanhatalo | Bayesian modeling with Gaussian processes using the MATLAB toolbox GPstuff[END_REF].

In practice, all implementations require the user to specify the mean and covariance functions of a Gaussian process prior under a parameterized form. Out of the various methods available to estimate the model parameters, we can safely say that the most popular approach is the maximum likelihood estimation (MLE) method. However, a simple numerical experiment consisting in interpolating a function (see Table 1), as Table 1: Inconsistencies in the results across different Python packages. The results were obtained by fitting a GP model, with constant mean and a Matérn kernel (ν = 5/2), to the Branin function, using the default settings for each package. We used 50 training points and 500 test points sampled from a uniform distribution on [-5, 10] × [0, 15]. The table reports the estimated values for the variance and length scale parameters of the kernel, the empirical root mean squared prediction error (ERMSPE) and the minimized negative log likelihood (NLL). The last row shows the improvement using the recommendations in this study. is usually done in Bayesian optimization, shows that different MLE implementations from different Python packages produce very dispersed numerical results when the default settings of each implementation are used. These significant differences were also noticed by [START_REF] Erickson | Comparison of Gaussian process modeling software[END_REF] but the causes and possible mitigation were not investigated. Note that each package uses its own default algorithm for the optimization of the likelihood: GPyTorch uses ADAM [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF], OpenTURNS uses a truncated Newton method [START_REF] Nash | Newton-type minimization via the Lanczos method[END_REF] and the others generally use L-BFGS-B [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF]. It turns out that none of the default results in Table 1 are really satisfactory compared to the result obtained using the recommendations in this study3 . Focusing on the case of GP interpolation (with Bayesian optimization as the main motivation), the first contribution of this article is to understand the origin of the inconsistencies across available implementations. The second contribution is to investigate simple but effective strategies for improving these implementations, using the wellestablished GPy package as a case study. We shall propose recommendations concerning several optimization settings: initialization and restart strategies, parameterization of the covariance, etc. By anticipation of our numerical results, the reader is invited to refer to Figure 1 and Table 2, which show that significant improvement in terms of estimated parameter values and prediction errors can be obtained over default settings using better optimization schemes.

LIBRARY

Even though this work targets a seemingly prosaic issue, and advocates somehow simple solutions, we feel that the contribution is nonetheless of significant value considering the widespread use of GP modeling. Indeed, a host of studies, particularly in the literature of Bayesian optimization, rely on off-the-shelf GP implementations: for their conclusions to be reliable and reproducible, robust implementations are critical.

The article is organized as follows. Section 2 provides a brief review of GP modeling and MLE. Section 3 describes some numerical aspects of the evaluation and optimization of the likelihood function, with a focus on GPy's implementation. Section 4 provides We remove one point at a time to obtain (a) the distribution of the differences of negative log-likelihood (NLL) values between the two setups; (b) the empirical CDFs of the prediction error at the removed points; (c) pairs of box-plots for the estimated range parameters (for each dimension, indexed from 1 to 8 on the x-axis, the box-plot for improved setup is on the left and the box-plot for default setup is on the right; horizontal red lines correspond to the estimated values using the whole data set without leave-one-out). Notice that the parameter distributions of the default setup are more spread out.

× R d → R. The objective is to predict Z(x) at a given location x ∈ R d , given a data set D = {(x i , z i ) ∈ R d × R, 1 ≤ i ≤ n},
where the observations z i s are assumed to be the outcome of an additive-noise model:

Z i = Z(x i ) + ε i , 1 ≤ i ≤ n.
In most applications, it is assumed that the ε i s are zero-mean Gaussian i.i.d. random variables with variance σ 2 ε ≥ 0, independent of Z. (In rarer cases, heteroscedasticity is assumed.) Knowing m and k, recall (see, e.g. [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] that the posterior distribution of Z is such that Z | Z 1 , . . . , Z n , m, k ∼ GP( Ẑn , k n ), where Ẑn and k n stand respectively for the posterior mean and covariance functions:

Ẑn (x) = m(x) + ∑ n i=1 w i (x; x n ) (z i -m(x i )) , k n (x, y) = k(x, y) -w(y; x n ) T K(x n , x) ,
where x n denotes observation points (x 1 , . . . , x n ) and the weights w i (x; x n ) are solutions of the linear system:

(K(x n , x n ) + σ 2 ε I n ) w(x; x n ) = K(x n , x) , (1)
with K(x n , x n ) the n × n covariance matrix with entries k(x i , x j ), I n the identity matrix of size n, and w(x; x n ) (resp. K(x n , x)) the column vector with entries w i (x;

x n ) (resp. k(x i , x)), 1 ≤ i ≤ n.
The posterior covariance at x, y ∈ R d may be written as

k n (x, y) = k(x, y) -w(y; x n ) T K(x n , x) . (2)
It is common practice to assume a zero mean function m = 0-a reasonable choice if the user has taken care to center data-but most GP implementations also provide an Table 3: Some kernel functions available in GPy. The Matérn kernel is recommended by [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF]. Γ denotes the gamma function, K ν is the modified Bessel function of the second kind.

KERNEL

r(h), h ∈ [0, +∞)

SQUARED EXPONENTIAL exp(-1 2 r 2 ) RATIONAL QUADRATIC (1 + r 2 ) -ν MAT ÉRN WITH PARAM. ν > 0 2 1-ν Γ (ν) √ 2νr ν K ν √ 2νr
option for setting a constant mean function m( • ) = µ ∈ R. In this article, we will include such a constant in our models, and treat it as an additional parameter to be estimated by MLE along with the others. (Alternatively, µ could be endowed with a Gaussian or improper-uniform prior, and then integrated out; see, e.g., [START_REF] O'hagan | Curve fitting and optimal design for prediction[END_REF].)

The covariance function, aka covariance kernel, models similarity between data points and reflects the user's prior belief about the function to be learned. Most GP implementations provide a couple of stationary covariance functions taken from the literature (e.g., [START_REF] Wendland | Scattered data approximation[END_REF][START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. The squared exponential, the rational quadratic or the Matérn covariance functions are popular choices (see Table 3). These covariance functions include a number of parameters: a variance parameter σ 2 > 0 corresponding to the variance of Z, and a set of range (or length scale) parameters ρ 1 , . . . , ρ d , such that

k(x, y) = σ 2 r(h) , (3)
with

h 2 = ∑ d i=1 (x [i] -y [i] ) 2 /ρ 2 i , where x [i]
and y [i] denote the elements of x and y. The function r : R → R in (3) is the stationary correlation function of Z. From now on, the vector of model parameters will be denoted by θ = (σ 2 , ρ 1 , . . . , ρ d , . . . , σ 2 ε ) T ∈ Θ ⊂ R p , and the corresponding covariance matrix K(x n , x n ) + σ 2 ε I n by K θ .

Maximum likelihood estimation

In this article, we focus on GP implementations where the parameters (θ , µ) ∈ Θ × R of the process Z are estimated by maximizing the likelihood L(Z n |θ , µ) of Z n = (Z 1 , . . . , Z n ) T , or equivalently, by minimizing the negative log-likelihood (NLL)

-log(L(Z n |θ , µ)) = 1 2 (Z n -µ1 n ) K -1 θ (Z n -µ1 n ) + 1 2 log|K θ | + constant. ( 4
)
This optimization is typically performed by gradient-based methods, although local maxima can be of significant concern as the likelihood is often non-convex. Computing the likelihood and its gradient with respect to (θ , µ) has a O(n 3 + dn 2 ) computational cost [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF][START_REF] Petit | Towards new cross-validation-based estimators for Gaussian process regression: efficient adjoint computation of gradients[END_REF].

Numerical noise

The evaluation of the NLL as well as its gradient is subject to numerical noise, which can prevent proper convergence of the optimization algorithms. Figure 2 shows a typical situation where the gradient-based optimization algorithm stops before converging to an actual minimum. In this section, we provide an analysis on the numerical noise on the NLL using the concept of local condition numbers. We also show that the popular solution of adding jitter cannot be considered as a fully satisfactory answer to the problem of numerical noise.

Numerical noise stems from both terms of the NLL, namely 1 2 Z n K -1 θ Z n and 1 2 log|K θ |. (For simplification, we assume µ = 0 in this section.)

First, recall that the condition number κ(K θ ) of K θ , defined as the ratio |λ max /λ min | of the largest eigenvalue to the smallest eigenvalue [START_REF] Press | Numerical recipes in C. The art of scientific computing[END_REF], is the key element for analyzing the numerical noise on K -1 θ Z n . In double-precision floating-point approximations of numbers, Z n is corrupted by an error ε whose magnitude is such that ε / Z n 10 -16 . Worst-case alignment of Z n and ε with the eigenvectors of K θ gives

K -1 θ ε K -1 θ Z n κ(K θ ) × 10 -16 , (5)
which shows how the numerical noise is amplified when K θ becomes ill-conditioned. The term log|K θ | is nonlinear in K θ , but observe, using the identity

d log|K θ |/dK θ = K -1 θ , that the differential of log| • | at K θ is given by H → Trace(K -1 θ H).
Thus, the induced operator norm with respect to the Frobenius norm • F is K -1 θ F . We can then apply results from [START_REF] Trefethen | Numerical Linear Algebra[END_REF] to get a local condition number of the mapping A → log|A| at K θ :

κ(log| • |, K θ ) lim ε→0 sup δ A F ≤ε log|K θ + δ A | -log|K θ | log|K θ | K θ F δ A F = ∑ n i=1 1 λ 2 i ∑ n i=1 λ 2 i |∑ n i=1 log(λ i )| (6) where λ 1 , • • • , λ n are the (positive) eigenvalues of K θ . Then, we have κ(K θ ) |∑ n i=1 log(λ i )| ≤ κ(log| • |, K θ ) ≤ nκ(K θ ) |∑ n i=1 log(λ i )| , (7)
which shows that numerical noise on log|K θ | is linked to the condition number of K θ .

The local condition number of the quadratic form 1 2 Z T n K -1 θ Z n as a function of Z n can also be computed analytically. Some straightforward calculations show that it is bounded by κ(K θ ).

(When the optimization algorithm stops in the example of Figure 2, we have κ(K θ ) 10 11 and κ(log| • |, K θ ) 10 9.5 . The empirical numerical fluctuations are measured as the residuals of a local second-order polynomial best fit, giving noise levels 10 -7 , 10 -8 and 10 -7.5 for K -1 θ Z n , 1 2 Z T n K -1 θ Z n and log|K θ | respectively. These values are consistent with the above first-order analysis.) Thus, when κ(K θ ) becomes large in the course of the optimization procedure, numerical noise on the likelihood and its gradient may trigger an early stopping of the optimization algorithm (supposedly when the algorithm is unable to find a proper direction of improvement). It is well-known that κ(K θ ) becomes large when σ 2 ε = 0 and one of the following conditions occurs: 1) data points are close, 2) the covariance is very smooth (as for instance when considering the squared exponential covariance), 3) when the range parameters ρ i are large. These conditions arise more often than not. Therefore, the problem of numerical noise in the evaluation of the likelihood and its gradient is a problem that should not be neglected in GP implementations.

The most classical approach to deal with ill-conditioned covariance matrices is to add a small positive number on the diagonal of the covariance matrix, called jitter, which is equivalent to assuming a small observation noise with variance σ 2 ε > 0. In GPy for instance, the strategy consists in always setting a minimal jitter of 10 -8 , which is automatically increased by an amount ranging from 10 -6 σ 2 to 10 -1 σ 2 whenever the Cholesky factorization of the covariance matrix fails (due to numerical non-positiveness). The smallest jitter making K θ numerically invertible is kept and an error is thrown if no jitter allows for successful factorization. However, note that large values for the jitter may yield smooth, non-interpolating approximations, with possible unintuitive and undesirable effects (see [START_REF] Andrianakis | The effect of the nugget on Gaussian process emulators of computer models[END_REF], and causing possible convergence problems in Bayesian optimization.

Table 4 illustrates the behaviour of GP interpolation when σ 2 ε is increased. It appears that finding a satisfying trade-off between good interpolation properties and low numerical noise level can be difficult. Table 4 also supports the connection in ( 5) and ( 7) between noise levels and κ(K θ ). In view of the results of Figure 1 based on the default settings of GPy and Table 4, we believe that adaptive jitter cannot be considered as a do-it-all solution.

Strategies for improving likelihood maximization

In this section we investigate simple but hopefully efficient levers / strategies to improve available implementations of MLE for GP interpolation, beyond the control of the numerical noise on the likelihood using jitter. We mainly focus on 1) initialization methods for the optimization procedure, 2) stopping criteria, 3) the effect of "restart" strategies and 4) the effect of the parameterization of the covariance.

Table 4: Influence of the jitter on the GP model (same setting as in Figure 2). The table reports the condition numbers κ(K θ ) and κ(log| • |, K θ ), and the impact on the relative empirical standard deviations δ quad and δ logdet of the numerical noise on Z T n K -1 θ Z n and log|K θ | respectively (measured using second-order polynomial regressions). As σ ε increases, δ quad and δ logdet decrease but the interpolation error SSR/SST = 1 n ∑ n j=1 (Z j -Ẑn (x j )) 2 /std(Z 1 , ..., Z n ) and the NLL increase. Reducing numerical noise while keeping good interpolation properties requires careful attention in practice.

σ 2 ε / σ 2 0.0 10 -8 10 -6 10 -4 10 -2 κ(K θ ) 10 11 10 9 10 7.5 10 5.5 10 3.5 κ(log| • |, K θ ) 10 9.5 10 8.5 10 6.5 10 4.5 10 2.5 δ quad 10 -8 (= 10 11-19 ) 10 -9.5 (= 10 9-18.5 ) 10 -10.5 (= 10 7.5-18 ) 10 -12 (= 10 5.5-17.5 ) 10 -14 (= 10 3.5-17.5 ) δ logdet 10 -7.5 (= 10 9.5-17 ) 10 -9 (= 10 8.5-17.5 ) 10 -11 (= 10 6.5-17.5 ) 10 -13.5 (= 10 4.5-18 ) 10 -15.5 (= 10 2.

Initialization strategies

Most GP implementations use a gradient-based local optimization algorithm to maximize the likelihood that requires the specification of starting/initial values for the parameters.

In the following, we consider different initialization strategies.

Moment-based initialization. A first strategy consists in setting the parameters using empirical moments of the data. More precisely, assuming a constant mean m = µ, and a stationary covariance k with variance σ 2 and range parameters ρ 1 , . . . , ρ d , set

µ init = mean (Z 1 , . . . , Z n ), (8)
σ 2 init = var (Z 1 , . . . , Z n ), (9)
ρ k, init = std (x 1, [k] , . . . , x n, [k] ), k = 1, . . . , d, (10)
where mean, var and std stand for the empirical mean, variance and standard deviation, and x i, [k] denotes the kth coordinate of x i ∈ R d . The rationale behind (10) (following, e.g., [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] is that the range parameters can be thought of as the distance one has to move in the input space for the function value to change significantly and we assume, a priori, that this distance is linked to the dispersion of data points.

In GPy for instance, the default initialization consists in setting µ = 0, σ 2 = 1 and ρ k = 1 for all k. This is equivalent to the moment-based initialization scheme when the data (both inputs and outputs) are centered and standardized. The practice of standardizing the input domain into a unit length hypercube has been proposed (see, e.g., [START_REF] Snoek | Practical Bayesian optimization of machine learning algorithms[END_REF] to deal with numerical issues that arise due to large length scale values.

Profiled initialization. Assume the range parameters ρ 1 , . . . , ρ d (and more generally, all parameters different from σ 2 , σ 2 ε and µ) are fixed, and set σ 2 ε = ασ 2 , with a prescribed multiplicative factor α ≥ 0. In this case, the NLL can be optimized analytically w.r.t. µ and σ 2 . Optimal values turn out to be the generalized least squares solutions

µ GLS = (1 T n K -1 θ 1 n ) -1 1 T n K -1 θ Z n , (11)
σ 2 GLS = 1 n (Z n -µ GLS 1 n ) T K -1 θ (Z n -µ GLS 1 n ) , (12)
where θ = (σ 2 , ρ 1 , . . . , ρ d , . . . , σ 2 ε ) T ∈ Θ , with σ 2 = 1 and σ 2 ε = α. Under the profiled initialization scheme, ρ 1 , . . . , ρ d are set using (10), α is prescribed according to user's preference, and µ and σ 2 are initialized using ( 11) and ( 12).

Grid-search initialization. Grid-search initialization is a profiled initialization with the addition of a grid-search optimization for the range parameters.

Define a nominal range vector ρ 0 such that

ρ 0,[k] = √ d max 1≤i≤n x i,[k] -min 1≤i≤n x i,[k] , 1 ≤ k ≤ d.
Then, define a one-dimensional grid of size L (e.g., L = 5) by taking range vectors proportional to ρ 0 : {α 1 ρ 0 , . . . , α L ρ 0 }, where the α i s range, in logarithmic scale, from a "small" value (e.g., α 1 = 1/50) to a "large" value (e.g., α L = 2). For each point of the grid, the likelihood is optimized with respect to µ and σ 2 using ( 11) and ( 12). The range vector with the best likelihood value is selected. (Note that this initialization procedure is the default initialization procedure in the Matlab/GNU Octave toolbox STK.)

Stopping condition

Most GP implementations rely on well-tested gradient-based optimization algorithms.

For instance, a popular choice in Python implementations is to use the limited-memory BFGS algorithm with box constraints (L-BFGS-B; see [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] of the SciPy ecosystem. (Other popular optimization algorithms include the ordinary BFGS, truncated Newton constrained, SQP, etc.; see, e.g., [START_REF] Nocedal | Numerical Optimization[END_REF].) The L-BFGS-B algorithm, which belongs to the class of quasi-Newton algorithms, uses limited-memory Hessian approximations and shows good performance on non-smooth functions [START_REF] Curtis | A quasi-Newton algorithm for nonconvex, nonsmooth optimization with global convergence guarantees[END_REF]. Regardless of which optimization algorithm is chosen, the user usually has the possibility to tune the behavior of the optimizer, and in particular to set the stopping condition. Generally, the stopping condition is met when a maximum number of iterations is reached or when a norm on the steps and/or the gradient become smaller than a threshold.

By increasing the strictness of the stopping condition during the optimization of the likelihood, one would expect better parameter estimations, provided the numerical noise on the likelihood does not interfere too much.

Restart and multi-start strategies

Due to numerical noise and possible non-convexity of the likelihood with respect to the parameters, gradient-based optimization algorithms may stall far from the global optimum. A common approach to circumvent the issue is to carry out several optimization runs with different initialization points. Two simple strategies can be compared.

Table 5: Two popular reparameterization mappings τ, as implemented, for example, in GPy and STK respectively. For invsoftplus, notice parameter s > 0, which is introduced when input standardization is considered (see Section 5).

REPARAM. METHOD

τ : R + → R τ -1 : R → R + INVSOFTPLUS(s) log(exp(θ /s) -1) s log(exp(θ ) + 1) LOG log(θ ) exp(θ )
Restart. In view of Figure 2, a first simple strategy is to restart the optimization algorithm to clear its memory (Hessian approximation, step sizes. . . ), hopefully allowing it to escape a possibly problematic location using the last best parameters as initial values for the next optimization run. The optimization can be restarted a number of times, until a budget N opt of restarts is spent or the best value for the likelihood does not improve.

Multi-start. Given an initialization point (θ init , µ init ) ∈ Θ × R, a multi-start strategy consists in running N opt > 1 optimizations with different initialization points corresponding to perturbations of the initial point (θ init , µ init ). In practice, we suggest the following rule for building the perturbations: first, move the range parameters around (ρ 1, init , . . . , ρ d, init ) T (refer to Section 5 for an implementation); then, propagate the perturbations on µ and σ 2 using ( 11) and ( 12). The parameter with the best likelihood value over all optimization runs is selected.

Parameterization of the covariance function

The parameters of the covariance functions are generally positive real numbers (σ 2 , ρ 1 , ρ 2 . . .) and are related to scaling effects that act "multiplicatively" on the predictive distributions. Most GP implementations introduce a reparameterization using a monotonic one-to-one mapping τ : R + → R, acting component-wise on the positive parameters of θ , resulting in a mapping τ : Θ → Θ . Thus, for carrying out MLE, the actual criterion J that is optimized in most implementations may then be written as

J : θ ∈ Θ → -log(L(Z n |τ -1 (θ ), c)) . (13)
Table 5 lists two popular reparameterization mappings τ.

The effect of reparameterization is to "reshape" the likelihood. Typical likelihood profiles using the log and the so-called invsoftplus reparameterizations are shown on Figure 3. Notice that the NLL may be almost flat in some regions depending on the reparameterization. Changing the shape of the optimization criterion, combined with numerical noise, may or may not facilitate the convergence of the optimization.

5 Numerical study

Methodology

The main metric used in this numerical study is based on empirical cumulative distributions (ECDFs) of differences on NLL values. More precisely, consider N + 1 optimization schemes S 0 , S 1 , . . . , S N , where S 0 stands for a "brute-force" optimization scheme based on a very large number of multi-starts, which is assumed to provide a robust MLE, and S 1 , . . . , S N are optimization schemes to be compared. Each optimization scheme is run on M data sets D j , 1 ≤ j ≤ M, and we denote by e i, j the difference

e i, j = NLL i, j -NLL 0, j , 1 ≤ i ≤ N, 1 ≤ j ≤ M,
where NLL i, j the NLL value obtained by optimization scheme S i on data set D j .

A good scheme S i should concentrate the empirical distribution of the sample E i = {e i, j , j = 1, . . . , M} around zero-in other words, the ECDF is close to the ideal CDF e → 1 [0,∞[ (e). Using ECDF also provides a convenient way to compare performances: a strategy with a "steeper" ECDF, or larger area under the ECDF, is better.

Optimization schemes

All experiments are performed using GPy version 1.9.9, with the default L-BFGS-B algorithm. We use a common setup and vary the configurations of the optimization levers as detailed below.

Common setup. All experiments use an estimated constant mean-function, an anisotropic Matérn covariance function with regularity ν = 5/2, and we assume no observation noise (the adaptive jitter of GPy ranging from 10 -6 σ 2 to 10 2 σ 2 is used, however).

Initialization schemes. Three initialization procedures from Section 4.1 are considered.

Stopping criteria. We consider two settings for the stopping condition of the L-BFGS-B algorithm, called soft (the default setting: maxiter= 1000, factr=10 7 , pgtol10 -5 ) and strict (maxiter= 1000, factr=10, pgtol= 10 -20 ).

Restart and multi-start. The two strategies of Section 4.3 are implemented using a log reparameterization and initialization points (θ init , µ init ) determined using a grid-search strategy. For the multi-start strategy the initial range parameters are perturbed according to the rule ρ ← ρ init • 10 η where η is drawn from a N (0, σ 2 η ) distribution. We take σ η = log 10 (5)/1.96 (≈ 0.35), to ensure that about 0.95 of the distribution of ρ is in the interval

[1/5 • ρ init , 5 • ρ init ].
Reparameterization. We study the log reparameterization and two variants of the invsoftplus. The first version called no-input-standardization simply corresponds to taking s = 1 for each range parameter. The second version called input-standardization consists in scaling the inputs to a unit standard deviation on each dimension (by taking the corresponding value for s).

Data sets

The data sets are generated from six well-known test functions in the literature of Bayesian optimization: the Branin function (d = 2; see, e.g. [START_REF] Surjanovic | Virtual library of simulation experiments: Test functions and datasets[END_REF], the Borehole function (d = 8; see, e.g. [START_REF] Worley | Deterministic uncertainty analysis[END_REF], the Welded Beam Design function (d = 4; see [START_REF] Chafekar | Constrained multi-objective optimization using steady state genetic algorithms[END_REF], the g10 function (d = 8; see Ahmed, 2004, p. 128), along with two modified versions, g10mod and g10modmod (see [START_REF] Feliot | Une approche bayésienne pour l'optimisation multi-objectif sous contrainte[END_REF].

Each function is evaluated on Latin hypercube samples with a multi-dimensional uniformity criterion (LHS-MDU; [START_REF] Deutsch | Latin hypercube sampling with multidimensional uniformity[END_REF], with varying sample size n ∈ {3d, 5d, 10d, 20d}, resulting in a total of 6 × 4 = 24 data sets.

Results and findings

Figure 4 shows the effect of reparameterization and the initialization method. Observe that the log reparameterization performs significantly better than the invsoftplus reparameterizations. For the log reparameterization, observe that the grid-search strategy brings a moderate but not negligible gain with respect to the two other initialization strategies, which behave similarly.

Next, we study the effect of the different restart strategies and the stopping conditions, on the case of the log reparameterization and grid-search initialization. The metric used for the comparison is the area under the ECDFs of the differences of NLLs, computed by integrating the ECDF between 0 and NLL max = 100. Thus, a perfect optimization strategy would achieve an area under the ECDF equal to 100. Since the multi-start strategy is stochastic, results are averaged over 50 repetitions of the optimization procedures (for each N opt value, the optimization strategy is repeated 50 times). The areas are plotted against the computational run time. Run times are averaged over the repetitions in the case of the multi-start strategy.

Figure 5 shows that the soft stopping condition seems uniformly better. The restart strategy yields small improvements using moderate computational overhead. The multistart strategy is able to achieve the best results at the price of higher computational costs.

Conclusions and recommendations

Our numerical study has shown that the parameterization of the covariance function has the most significant impact on the accuracy of MLE in GPy. Using restart / multi-start strategies is also very beneficial to mitigate the effect of the numerical noise on the likelihood. The two other levers have second-order but nonetheless measurable influence. These observations make it possible to devise a recommended combination of improvement levers-for GPy at least, but hopefully transferable to other software packages as well. When computation time matters, an improved optimization procedure for MLE consists in choosing the combination of a log reparameterization, with a gridsearch initialization, the soft (GPy's default) stopping condition, and a small number, say N opt = 5, of restarts.

Figure 1 and Table 2 are based on the above optimization procedure, which results in significantly better likelihood values and smaller prediction errors. The multi-start strategy can be used when accurate results are sought.

Several topics could be investigated in the future: the optimization of alternative selection criteria, in particular criteria based on leave-one-out procedures, the case of regression, the problem of parameterization in relation to that of the identifiability of the parameters (see, e.g., [START_REF] Anderes | On the consistent separation of scale and variance for Gaussian random fields[END_REF].

As a conclusion, our recommendations are not intended to be universal, but will hopefully encourage researchers and users to develop and use more reliable and more robust GP implementations, in Bayesian optimization or elsewhere.

Figure 1 :

 1 Figure1: Improved (cf. Section 6) vs default setups in GPy on the Borehole function with n = 20d = 160 random training points. We remove one point at a time to obtain (a) the distribution of the differences of negative log-likelihood (NLL) values between the two setups; (b) the empirical CDFs of the prediction error at the removed points; (c) pairs of box-plots for the estimated range parameters (for each dimension, indexed from 1 to 8 on the x-axis, the box-plot for improved setup is on the left and the box-plot for default setup is on the right; horizontal red lines correspond to the estimated values using the whole data set without leave-one-out). Notice that the parameter distributions of the default setup are more spread out.

Figure 2 :

 2 Figure2: Noisy NLL profile along a particular direction in the parameter space, with a best linear fit (orange line). This example was obtained with GPy while estimating the parameters of a Matérn 5/2 covariance, using 20 data points sampled from a Branin function, and setting σ 2 ε = 0. The red vertical line indicates the location where the optimization of the likelihood stalled.

Figure 3 :

 3 Figure3: Profiles of the NLL along a linear path t through the profiled initialization point (at zero, blue vertical line) and the optimum (at one, black vertical line). Orange (resp. blue) line corresponds to the log (resp. invsoftplus) reparameterization.

Figure 4 :Figure 5 :

 45 Figure 4: Initialization and reparameterization methods. (a) ECDFs corresponding to the best initialization method for each of the three reparameterizations-red line: log reparam. with gridsearch init.; green line: invsoftplus with input-standardization reparam. and grid-search init; blue line: invsoftplus with no-input-standardization reparam. and moment-based init. (b) ECDFs for different initialization methods for the log reparameterization.

Table 2 :

 2 Improved (cf. Section 6) vs default setups in GPy for the interpolation of the Borehole

	function (input space dimension is d = 8) with n ∈ {3d, 5d} random data points (see Section 5.3
	for details). The experiment is repeated 50 times. The columns report the leave-one-out mean
	squared error (LOO-MSE) values (empirical mean over the repetitions, together with the standard
	deviation and the average proportion of the LOO-MSE to the total standard deviation of the data
	in parentheses).
	METHOD	n = 3d	n = 5d
	DEFAULT	17.559 (4.512, 0.387)	10.749 (2.862, 0.229)
	IMPROVED	3.949 (1.447, 0.087)	1.577 (0.611, 0.034)
	an analysis of factors influencing the accuracy of numerical MLE procedures. Finally,
	Section 5 assesses the effectiveness of our solutions through numerical experiments and
	Section 6 concludes the article.
	2 Background
	2.1 Gaussian processes
	Let Z ∼ GP(m, k) be a Gaussian process indexed by R d , d ≥ 1, specified by a mean
	function m : R

d → R and a covariance function k : R d

Code available at https://github.com/saferGPMLE